HCC: US vs. UK

library(ggplot2)
library(survminer)

## Loading required package: ggpubr

library(survival)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ tibble  3.1.0     ✓ purrr   0.3.4
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

#library(data.table) 
library(forestmodel) 
#library(Hmisc) 
library(sjPlot)

## Registered S3 methods overwritten by 'lme4':
##   method                          from
##   cooks.distance.influence.merMod car 
##   influence.merMod                car 
##   dfbeta.influence.merMod         car 
##   dfbetas.influence.merMod        car

library(stargazer)

## 
## Please cite as:

##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.

##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer

library(sjmisc)

## 
## Attaching package: 'sjmisc'

## The following object is masked from 'package:purrr':
## 
##     is_empty

## The following object is masked from 'package:tidyr':
## 
##     replace_na

## The following object is masked from 'package:tibble':
## 
##     add_case

library(arsenal)

## 
## Attaching package: 'arsenal'

## The following object is masked from 'package:sjmisc':
## 
##     %nin%

library(gtsummary)
library(expss)

## 
## Use 'expss_output_viewer()' to display tables in the RStudio Viewer.
##  To return to the console output, use 'expss_output_default()'.

## 
## Attaching package: 'expss'

## The following objects are masked from 'package:gtsummary':
## 
##     contains, vars

## The following objects are masked from 'package:sjmisc':
## 
##     add_columns, add_rows, rec

## The following objects are masked from 'package:stringr':
## 
##     fixed, regex

## The following objects are masked from 'package:purrr':
## 
##     keep, modify, modify_if, transpose, when

## The following objects are masked from 'package:tidyr':
## 
##     contains, nest

## The following objects are masked from 'package:dplyr':
## 
##     between, compute, contains, first, last, na_if, recode, vars

## The following object is masked from 'package:ggpubr':
## 
##     compare_means

## The following object is masked from 'package:ggplot2':
## 
##     vars

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following object is masked from 'package:arsenal':
## 
##     is.Date

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(ggsignif)
library(ggsci)
library(Greg) #to use timesplitter

## Loading required package: forestplot

## Loading required package: grid

## Loading required package: magrittr

## 
## Attaching package: 'magrittr'

## The following objects are masked from 'package:expss':
## 
##     and, equals, not, or

## The following object is masked from 'package:arsenal':
## 
##     set_attr

## The following object is masked from 'package:purrr':
## 
##     set_names

## The following object is masked from 'package:tidyr':
## 
##     extract

## Loading required package: checkmate

## Loading required package: Gmisc

## Loading required package: Rcpp

## Loading required package: htmlTable

## Registered S3 methods overwritten by 'Hmisc':
##   method                 from 
##   [.labelled             expss
##   print.labelled         expss
##   as.data.frame.labelled expss

library(magrittr)
options(scipen=999)
#Random effects
#https://stats.idre.ucla.edu/r/dae/mixed-effects-cox-regression/

library("coxme")

## Loading required package: bdsmatrix

## 
## Attaching package: 'bdsmatrix'

## The following object is masked from 'package:base':
## 
##     backsolve

Read in data

bound2 <- read_csv("/Users/Ivanics/Desktop/GitHub/International-LDLT/Datafile/merged.csv", guess_max = 300000)

## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   .default = col_double(),
##   DTYPE = col_character(),
##   DONCOD = col_character(),
##   DCMV = col_character(),
##   DSEX = col_character(),
##   BLD_GP_MATCH = col_character(),
##   GRAFT_TYPE = col_character(),
##   RSEX = col_character(),
##   RETHNIC = col_character(),
##   TRANSPLANT_UNIT = col_character(),
##   RREN_SUP = col_character(),
##   RVENT = col_character(),
##   RAB_SURGERY = col_character(),
##   RLIFE = col_character(),
##   RASCITES = col_character(),
##   RENCEPH = col_character(),
##   RBG = col_character(),
##   RANTI_HCV = col_character(),
##   COUNTRY = col_character(),
##   HCC_combined = col_character(),
##   UKT_PLDGRP = col_character()
##   # ... with 5 more columns
## )
## ℹ Use `spec()` for the full column specifications.

bound2 <- bound2 %>% filter(COUNTRY == "US" | COUNTRY == "UK" | COUNTRY == "CAN")

bound2 <- bound2 %>% filter(TX_YR <= 2018)

bound2 <- bound2 %>%
  mutate(PSURV_years = PSURV/365.25) %>%
  mutate(GSURV_years = GSURV/365.25)

bound2$COUNTRY <- factor(bound2$COUNTRY, ordered = FALSE)
bound2$COUNTRY <- relevel(bound2$COUNTRY, ref="US")
bound2$RVENT <- factor(bound2$RVENT, ordered = FALSE)
bound2$RVENT <- relevel(bound2$RVENT, ref = "Not ventilated")
bound2$RETHNIC <- factor(bound2$RETHNIC, ordered = FALSE)
bound2$RETHNIC <- relevel(bound2$RETHNIC, ref = "White")
bound2$DCMV <- factor(bound2$DCMV, ordered = FALSE)
bound2$DCMV <- relevel(bound2$DCMV, ref = "Negative")
bound2$HCC <- factor(bound2$HCC)
bound2$HCV <- factor(bound2$HCV)
bound2$ALD <- factor(bound2$ALD)
bound2$NASH <- factor(bound2$NASH)
bound2$PSC <- factor(bound2$PSC)
bound2$PBC <- factor(bound2$PBC)
bound2$MET <- factor(bound2$MET)
bound2$AID <- factor(bound2$AID)
bound2$TRANSPLANT_UNIT <- factor(bound2$TRANSPLANT_UNIT)

bound2 <- bound2 %>% mutate(LDLT = case_when(
  DTYPE == "LDLT" & GRAFT_TYPE == "Segmental" ~ 1,
  TRUE ~ 0
)) %>% mutate(LDLT = factor(LDLT))

bound2 <- bound2 %>% mutate(MELD_calculated = 3.78*log(RBILIRUBIN) + 11.2*log(RINR) + 9.57*log(RCREAT) + 6.43)

Recoding graft survival

#bound2 <- bound2 %>% mutate(
#  GCENSnew = case_when(
#    PCENS == 1 & PSURV_years <= GSURV_years ~ 1,
#    TRUE ~ 0
#  )
#)

#bound2 <- bound2 %>% mutate(
#  GSURV_yearsnew = case_when(
#    PCENS == 1 & PSURV_years <= GSURV_years ~ PSURV_years,
#    TRUE ~ GSURV_years
#  )
#)

Recoding patient survival

#1year PSURV
bound2 <- bound2 %>% mutate(PSURV_1year =
    case_when(
      PSURV_years >= 1.0000001 ~ 1,
      PSURV_years <= 1 ~ PSURV_years
    ))

#1-year PCENS
bound2 <- bound2 %>% mutate(PCENS_1year =
                              case_when(
                                PSURV_years <= 1 ~ PCENS,
                                PSURV_years > 1 ~ 0
                              ))

#3year PSURV
bound2 <- bound2 %>% mutate(PSURV_3year =
    case_when(
      PSURV_years >= 3.000001 ~ 3,
      PSURV_years <= 3 ~ PSURV_years
    ))

#3-year PCENS
bound2 <- bound2 %>% mutate(PCENS_3year =
                              case_when(
                                PSURV_years <= 3 ~ PCENS,
                                PSURV_years > 3 ~ 0
                              ))

#5year PSURV
bound2 <- bound2 %>% mutate(PSURV_5year =
    case_when(
      PSURV_years >= 5.000001 ~ 5,
      PSURV_years <= 5 ~ PSURV_years
    ))

#5-year PCENS
bound2 <- bound2 %>% mutate(PCENS_5year =
                              case_when(
                                PSURV_years <= 5 ~ PCENS,
                                PSURV_years > 5 ~ 0
                              ))

#10year PSURV
bound2 <- bound2 %>% mutate(PSURV_10year =
    case_when(
      PSURV_years >= 10.000001 ~ 10,
      PSURV_years <= 10 ~ PSURV_years
    ))

#5-year PCENS
bound2 <- bound2 %>% mutate(PCENS_10year =
                              case_when(
                                PSURV_years <= 10 ~ PCENS,
                                PSURV_years > 10 ~ 0
                              ))

Recoding graft survival

#1year GSURV
bound2 <- bound2 %>% mutate(GSURV_1year =
    case_when(
      GSURV_years >= 1.0000001 ~ 1,
      GSURV_years <= 1 ~ GSURV_years
    ))

#1-year GCENS
bound2 <- bound2 %>% mutate(GCENS_1year =
                              case_when(
                                GSURV_years <= 1 ~ GCENS,
                                GSURV_years > 1 ~ 0
                              ))

#3year GSURV
bound2 <- bound2 %>% mutate(GSURV_3year =
    case_when(
      GSURV_years >= 3.000001 ~ 3,
      GSURV_years <= 3 ~ GSURV_years
    ))

#3-year GCENS
bound2 <- bound2 %>% mutate(GCENS_3year =
                              case_when(
                                GSURV_years <= 3 ~ GCENS,
                                GSURV_years > 3 ~ 0
                              ))

#5year GSURV
bound2 <- bound2 %>% mutate(GSURV_5year =
    case_when(
      GSURV_years >= 5.000001 ~ 5,
      GSURV_years <= 5 ~ GSURV_years
    ))

#5-year GCENS
bound2 <- bound2 %>% mutate(GCENS_5year =
                              case_when(
                                GSURV_years <= 5 ~ GCENS,
                                GSURV_years > 5 ~ 0
                              ))

#10year GSURV
bound2 <- bound2 %>% mutate(GSURV_10year =
    case_when(
      GSURV_years >= 10.000001 ~ 10,
      GSURV_years <= 10 ~ GSURV_years
    ))

#5-year GCENS
bound2 <- bound2 %>% mutate(GCENS_10year =
                              case_when(
                                GSURV_years <= 10 ~ GCENS,
                                GSURV_years > 10 ~ 0
                              ))

bound3 <- bound2
bound3 <- bound3 %>% filter(!(DTYPE == "LDLT" & GRAFT_TYPE == "Whole"))
bound2 <- bound2 %>% filter(DTYPE == "LDLT" & GRAFT_TYPE == "Segmental")

#creatining a variable for center count (number of appearances of a center in LT pts)
bound2 <- bound2 %>% 
  group_by(TRANSPLANT_UNIT, TX_YR) %>%
  mutate(count = n()) %>%
  ungroup

summary(bound2$count)

Min. 1st Qu. Median Mean 3rd Qu. Max. 1.00 7.00 13.00 17.12 22.00 58.00

#creating a variable for center experience according to percentile of count
bound2 <- bound2 %>%
  mutate(Centerpercentile =
           case_when(
             count <7 ~ 1,
             count >=7 & count < 22 ~2,
             count >=22 ~3,
           )) %>%
  mutate(Centerpercentile = factor(Centerpercentile, labels = c("<25%ile", "25-75%ile", ">75%ile")))

Crosstabulating

#ftable(xtabs(cbind(PCENS,GCENS)~ COUNTRY, data=bound2))
#ftable(xtabs(cbind(PCENS,GCENSnew)~ COUNTRY, data=bound2))

bound3 <- bound3 %>% mutate(DTYPE_DDLTandLDLT = case_when(
  DTYPE == "DBD" ~ 0,
  DTYPE == "DCD" ~ 0,
  DTYPE == "LDLT" ~ 1
))  %>% mutate(DTYPE_DDLTandLDLT = factor(DTYPE_DDLTandLDLT, labels = c("DDLT", "LDLT")))

bound3_2008 <- bound3 %>% filter(TX_YR == 2008)
tab1 <- tableby(COUNTRY ~
              DTYPE_DDLTandLDLT,
                data=bound3_2008, test=TRUE, total=TRUE, 
                numeric.stats=c("medianq1q3"), numeric.test="kwt", cat.test="chisq")
summary(tab1, title='Graft types 2008', pfootnote=TRUE, digits = 2)

Graft types 2008
	US (N=4920)	CAN (N=284)	UK (N=473)	Total (N=5677)	p value
DTYPE_DDLTandLDLT					< 0.001¹
DDLT	4758 (96.7%)	227 (79.9%)	465 (98.3%)	5450 (96.0%)
LDLT	162 (3.3%)	57 (20.1%)	8 (1.7%)	227 (4.0%)

Pearson’s Chi-squared test

bound3_2018 <- bound3 %>% filter(TX_YR == 2018)
tab1 <- tableby(COUNTRY ~
              DTYPE_DDLTandLDLT,
                data=bound3_2018, test=TRUE, total=TRUE, 
                numeric.stats=c("medianq1q3"), numeric.test="kwt", cat.test="chisq")
summary(tab1, title='Graft types 2018', pfootnote=TRUE, digits = 2)

Graft types 2018
	US (N=6655)	CAN (N=348)	UK (N=771)	Total (N=7774)	p value
DTYPE_DDLTandLDLT					< 0.001¹
DDLT	6327 (95.1%)	301 (86.5%)	768 (99.6%)	7396 (95.1%)
LDLT	328 (4.9%)	47 (13.5%)	3 (0.4%)	378 (4.9%)

Pearson’s Chi-squared test

tab1 <- tableby(COUNTRY ~
              CIT+
              DAGE+
              DBMI+
              DTYPE+
              DSEX+
              GRAFT_TYPE+
              DTYPE+
              DCMV,
              data=bound2, test=TRUE, total=TRUE, 
                numeric.stats=c("medianq1q3"), numeric.test="kwt", cat.test="chisq")
summary(tab1, title='Donor characteristics', pfootnote=TRUE, digits = 2)

Donor characteristics
	US (N=2433)	CAN (N=556)	UK (N=97)	Total (N=3086)	p value
CIT					0.003¹
Median (Q1, Q3)	84.90 (55.80, 120.00)	87.00 (57.00, 127.00)	110.00 (72.00, 162.00)	86.00 (57.00, 123.00)
DAGE					0.010¹
Median (Q1, Q3)	36.00 (28.00, 45.00)	35.00 (27.00, 46.00)	32.50 (25.00, 40.00)	36.00 (28.00, 45.00)
DBMI
Median (Q1, Q3)	25.94 (23.46, 28.58)	25.34 (23.09, 28.33)	NA	25.90 (23.45, 28.57)
DTYPE					< 0.001²
LDLT	2433 (100.0%)	556 (100.0%)	97 (100.0%)	3086 (100.0%)
DSEX					< 0.001³
Female	1280 (52.6%)	242 (43.5%)	35 (36.1%)	1557 (50.5%)
Male	1153 (47.4%)	314 (56.5%)	62 (63.9%)	1529 (49.5%)
GRAFT_TYPE					< 0.001²
Segmental	2433 (100.0%)	556 (100.0%)	97 (100.0%)	3086 (100.0%)
DCMV
N-Miss	2433	433	97	2963
Negative	0	71 (57.7%)	0	71 (57.7%)
Positive	0	52 (42.3%)	0	52 (42.3%)

Kruskal-Wallis rank sum test
Chi-squared test for given probabilities
Pearson’s Chi-squared test

bound3 <- bound3 %>% mutate(UHN = factor(UHN, labels = c("Not UHN", "UHN")))

tab1 <- tableby(UHN ~
              DTYPE,
              data=bound3, subset = COUNTRY == "CAN", test=TRUE, total=TRUE, 
                numeric.stats=c("medianq1q3"), numeric.test="kwt", cat.test="chisq")
summary(tab1, title='', pfootnote=TRUE, digits = 2)

	Not UHN (N=2040)	UHN (N=1586)	Total (N=3626)	p value
DTYPE				< 0.001¹
DBD	1715 (84.1%)	1055 (66.5%)	2770 (76.4%)
DCD	190 (9.3%)	110 (6.9%)	300 (8.3%)
LDLT	135 (6.6%)	421 (26.5%)	556 (15.3%)

Pearson’s Chi-squared test

tab1 <- tableby(COUNTRY ~ 
              RSEX+
              RETHNIC+
              MELD+
              MELD_calculated+
              RREN_SUP+
              RAGE+
              RLIFE+
              RVENT+
              WAITLIST_TIME+
              TX_YR+
              BMI+
              RAB_SURGERY+
              RCREAT+
              RBILIRUBIN+
              RALBUMIN+
              RINR+
              RANTI_HCV+
              RBG+
              RASCITES+
              RENCEPH+ 
              BLD_GP_MATCH+
              HCC +
              HCV +
                ALD +
                NASH +
                PSC +
                PBC +
                AID+
                PSURV_years +
                factor(UHN),
                data=bound2, test=TRUE, total=TRUE, 
                numeric.stats=c("medianq1q3"), numeric.test="kwt", cat.test="chisq")
summary(tab1, title='Recipient characteristics', pfootnote=TRUE, digits = 2)

Recipient characteristics
	US (N=2433)	CAN (N=556)	UK (N=97)	Total (N=3086)	p value
RSEX					0.177¹
Female	1082 (44.5%)	242 (43.5%)	52 (53.6%)	1376 (44.6%)
Male	1351 (55.5%)	314 (56.5%)	45 (46.4%)	1710 (55.4%)
RETHNIC					< 0.001¹
White	1989 (81.8%)	123 (22.1%)	51 (52.6%)	2163 (70.1%)
Black	80 (3.3%)	1 (0.2%)	1 (1.0%)	82 (2.7%)
Other	364 (15.0%)	432 (77.7%)	45 (46.4%)	841 (27.3%)
MELD					< 0.001²
Median (Q1, Q3)	15.00 (11.00, 19.00)	17.00 (13.00, 21.00)	15.00 (11.00, 18.00)	15.00 (11.00, 19.00)
MELD_calculated					< 0.001²
Median (Q1, Q3)	12.65 (8.57, 17.03)	14.81 (9.99, 18.73)	11.07 (7.77, 15.64)	12.72 (8.65, 17.17)
RREN_SUP					0.013¹
N-Miss	5	0	6	11
No pre-tx support	2412 (99.3%)	549 (98.7%)	88 (96.7%)	3049 (99.2%)
Pre-tx support	16 (0.7%)	7 (1.3%)	3 (3.3%)	26 (0.8%)
RAGE					< 0.001²
Median (Q1, Q3)	56.00 (46.00, 62.00)	54.00 (45.00, 61.00)	53.00 (35.00, 61.00)	55.00 (46.00, 62.00)
RLIFE
N-Miss	31	556	1	588
High	742 (30.9%)	0	12 (12.5%)	754 (30.2%)
Intermediate	1258 (52.4%)	0	36 (37.5%)	1294 (51.8%)
Low	402 (16.7%)	0	48 (50.0%)	450 (18.0%)
RVENT
N-Miss	0	556	1	557
Not ventilated	2421 (99.5%)	0	96 (100.0%)	2517 (99.5%)
Ventilated	12 (0.5%)	0	0 (0.0%)	12 (0.5%)
WAITLIST_TIME					< 0.001²
Median (Q1, Q3)	151.00 (79.00, 308.00)	122.00 (63.50, 230.00)	91.50 (7.00, 201.75)	145.00 (75.00, 290.00)
TX_YR					< 0.001²
Median (Q1, Q3)	2014.00 (2011.00, 2017.00)	2013.00 (2010.00, 2015.00)	2013.00 (2011.00, 2015.00)	2014.00 (2011.00, 2016.00)
BMI					< 0.001²
Median (Q1, Q3)	26.25 (23.25, 29.86)	25.67 (22.67, 29.05)	24.24 (21.91, 28.20)	26.05 (23.11, 29.68)
RAB_SURGERY
N-Miss	37	556	1	594
No	1195 (49.9%)	0	78 (81.2%)	1273 (51.1%)
Yes	1201 (50.1%)	0	18 (18.8%)	1219 (48.9%)
RCREAT					0.058²
Median (Q1, Q3)	0.85 (0.70, 1.10)	0.82 (0.70, 1.05)	0.77 (0.57, 0.97)	0.84 (0.69, 1.10)
RBILIRUBIN					< 0.001²
Median (Q1, Q3)	2.70 (1.47, 4.80)	3.71 (1.99, 7.27)	2.49 (1.35, 4.50)	2.80 (1.50, 5.00)
RALBUMIN
Median (Q1, Q3)	3.10 (2.70, 3.60)	NA	3.20 (2.68, 3.82)	3.10 (2.70, 3.60)
RINR					0.039²
Median (Q1, Q3)	1.40 (1.20, 1.60)	1.44 (1.20, 1.79)	1.40 (1.20, 1.90)	1.40 (1.20, 1.66)
RANTI_HCV					0.760¹
N-Miss	76	111	0	187
Negative	1794 (76.1%)	339 (76.2%)	77 (79.4%)	2210 (76.2%)
Positive	563 (23.9%)	106 (23.8%)	20 (20.6%)	689 (23.8%)
RBG					< 0.001¹
N-Miss	0	1	0	1
0	1112 (45.7%)	237 (42.7%)	47 (48.5%)	1396 (45.3%)
A	1041 (42.8%)	206 (37.1%)	21 (21.6%)	1268 (41.1%)
AB	34 (1.4%)	25 (4.5%)	2 (2.1%)	61 (2.0%)
B	246 (10.1%)	87 (15.7%)	27 (27.8%)	360 (11.7%)
RASCITES
N-Miss	0	556	0	556
Ascites	1572 (64.6%)	0	44 (45.4%)	1616 (63.9%)
No ascites	861 (35.4%)	0	53 (54.6%)	914 (36.1%)
RENCEPH
N-Miss	0	556	3	559
Encephalopathic	1244 (51.1%)	0	22 (23.4%)	1266 (50.1%)
Not encephalopathic	1189 (48.9%)	0	72 (76.6%)	1261 (49.9%)
BLD_GP_MATCH
N-Miss	0	556	2	558
Compatible	539 (22.2%)	0	21 (22.1%)	560 (22.2%)
Identical	1880 (77.3%)	0	73 (76.8%)	1953 (77.3%)
Incompatible	14 (0.6%)	0	1 (1.1%)	15 (0.6%)
HCC					0.575¹
0	2023 (83.1%)	454 (81.7%)	78 (80.4%)	2555 (82.8%)
1	410 (16.9%)	102 (18.3%)	19 (19.6%)	531 (17.2%)
HCV					< 0.001¹
0	2383 (97.9%)	488 (87.8%)	77 (79.4%)	2948 (95.5%)
1	50 (2.1%)	68 (12.2%)	20 (20.6%)	138 (4.5%)
ALD					0.315¹
0	2173 (89.3%)	485 (87.2%)	88 (90.7%)	2746 (89.0%)
1	260 (10.7%)	71 (12.8%)	9 (9.3%)	340 (11.0%)
NASH					< 0.001¹
0	2029 (83.4%)	517 (93.0%)	86 (88.7%)	2632 (85.3%)
1	404 (16.6%)	39 (7.0%)	11 (11.3%)	454 (14.7%)
PSC					0.293¹
0	2037 (83.7%)	452 (81.3%)	78 (80.4%)	2567 (83.2%)
1	396 (16.3%)	104 (18.7%)	19 (19.6%)	519 (16.8%)
PBC					0.218¹
0	2239 (92.0%)	500 (89.9%)	87 (89.7%)	2826 (91.6%)
1	194 (8.0%)	56 (10.1%)	10 (10.3%)	260 (8.4%)
AID					0.076¹
0	2199 (90.4%)	516 (92.8%)	84 (86.6%)	2799 (90.7%)
1	234 (9.6%)	40 (7.2%)	13 (13.4%)	287 (9.3%)
PSURV_years					< 0.001²
Median (Q1, Q3)	3.40 (1.21, 6.05)	4.89 (2.40, 7.60)	3.18 (1.16, 5.50)	3.81 (1.52, 6.63)
factor(UHN)
N-Miss	2433	0	97	2530
0	0	135 (24.3%)	0	135 (24.3%)
1	0	421 (75.7%)	0	421 (75.7%)

Pearson’s Chi-squared test
Kruskal-Wallis rank sum test

tab1 <- tableby(COUNTRY ~ 
              factor(TRANSPLANT_UNIT),
                data=bound2, test=TRUE, total=TRUE, 
                numeric.stats=c("medianq1q3"), numeric.test="kwt", cat.test="chisq")
summary(tab1, title='Recipient characteristics', pfootnote=TRUE, digits = 2)

Recipient characteristics
	US (N=2433)	CAN (N=556)	UK (N=97)	Total (N=3086)	p value
factor(TRANSPLANT_UNIT)					< 0.001¹
0	0 (0.0%)	135 (24.3%)	0 (0.0%)	135 (4.4%)
00124	141 (5.8%)	0 (0.0%)	0 (0.0%)	141 (4.6%)
00248	1 (0.0%)	0 (0.0%)	0 (0.0%)	1 (0.0%)
00279	184 (7.6%)	0 (0.0%)	0 (0.0%)	184 (6.0%)
00403	1 (0.0%)	0 (0.0%)	0 (0.0%)	1 (0.0%)
00651	27 (1.1%)	0 (0.0%)	0 (0.0%)	27 (0.9%)
01395	3 (0.1%)	0 (0.0%)	0 (0.0%)	3 (0.1%)
01860	44 (1.8%)	0 (0.0%)	0 (0.0%)	44 (1.4%)
02573	88 (3.6%)	0 (0.0%)	0 (0.0%)	88 (2.9%)
02666	44 (1.8%)	0 (0.0%)	0 (0.0%)	44 (1.4%)
03503	7 (0.3%)	0 (0.0%)	0 (0.0%)	7 (0.2%)
03906	1 (0.0%)	0 (0.0%)	0 (0.0%)	1 (0.0%)
04464	5 (0.2%)	0 (0.0%)	0 (0.0%)	5 (0.2%)
05332	12 (0.5%)	0 (0.0%)	0 (0.0%)	12 (0.4%)
05549	13 (0.5%)	0 (0.0%)	0 (0.0%)	13 (0.4%)
05704	120 (4.9%)	0 (0.0%)	0 (0.0%)	120 (3.9%)
06107	9 (0.4%)	0 (0.0%)	0 (0.0%)	9 (0.3%)
06200	139 (5.7%)	0 (0.0%)	0 (0.0%)	139 (4.5%)
06603	19 (0.8%)	0 (0.0%)	0 (0.0%)	19 (0.6%)
07223	66 (2.7%)	0 (0.0%)	0 (0.0%)	66 (2.1%)
07347	74 (3.0%)	0 (0.0%)	0 (0.0%)	74 (2.4%)
07471	92 (3.8%)	0 (0.0%)	0 (0.0%)	92 (3.0%)
07905	17 (0.7%)	0 (0.0%)	0 (0.0%)	17 (0.6%)
08277	1 (0.0%)	0 (0.0%)	0 (0.0%)	1 (0.0%)
08587	1 (0.0%)	0 (0.0%)	0 (0.0%)	1 (0.0%)
09114	9 (0.4%)	0 (0.0%)	0 (0.0%)	9 (0.3%)
09362	47 (1.9%)	0 (0.0%)	0 (0.0%)	47 (1.5%)
1	0 (0.0%)	421 (75.7%)	0 (0.0%)	421 (13.6%)
10044	10 (0.4%)	0 (0.0%)	0 (0.0%)	10 (0.3%)
11067	83 (3.4%)	0 (0.0%)	0 (0.0%)	83 (2.7%)
11129	29 (1.2%)	0 (0.0%)	0 (0.0%)	29 (0.9%)
11191	93 (3.8%)	0 (0.0%)	0 (0.0%)	93 (3.0%)
11470	75 (3.1%)	0 (0.0%)	0 (0.0%)	75 (2.4%)
11532	1 (0.0%)	0 (0.0%)	0 (0.0%)	1 (0.0%)
11749	8 (0.3%)	0 (0.0%)	0 (0.0%)	8 (0.3%)
12834	3 (0.1%)	0 (0.0%)	0 (0.0%)	3 (0.1%)
13237	1 (0.0%)	0 (0.0%)	0 (0.0%)	1 (0.0%)
13485	1 (0.0%)	0 (0.0%)	0 (0.0%)	1 (0.0%)
13919	92 (3.8%)	0 (0.0%)	0 (0.0%)	92 (3.0%)
14477	39 (1.6%)	0 (0.0%)	0 (0.0%)	39 (1.3%)
15283	7 (0.3%)	0 (0.0%)	0 (0.0%)	7 (0.2%)
15438	6 (0.2%)	0 (0.0%)	0 (0.0%)	6 (0.2%)
15593	4 (0.2%)	0 (0.0%)	0 (0.0%)	4 (0.1%)
16523	1 (0.0%)	0 (0.0%)	0 (0.0%)	1 (0.0%)
16616	4 (0.2%)	0 (0.0%)	0 (0.0%)	4 (0.1%)
16864	8 (0.3%)	0 (0.0%)	0 (0.0%)	8 (0.3%)
17081	60 (2.5%)	0 (0.0%)	0 (0.0%)	60 (1.9%)
17918	21 (0.9%)	0 (0.0%)	0 (0.0%)	21 (0.7%)
18755	136 (5.6%)	0 (0.0%)	0 (0.0%)	136 (4.4%)
19034	8 (0.3%)	0 (0.0%)	0 (0.0%)	8 (0.3%)
20677	13 (0.5%)	0 (0.0%)	0 (0.0%)	13 (0.4%)
21080	182 (7.5%)	0 (0.0%)	0 (0.0%)	182 (5.9%)
21514	6 (0.2%)	0 (0.0%)	0 (0.0%)	6 (0.2%)
22692	50 (2.1%)	0 (0.0%)	0 (0.0%)	50 (1.6%)
24304	232 (9.5%)	0 (0.0%)	0 (0.0%)	232 (7.5%)
24335	9 (0.4%)	0 (0.0%)	0 (0.0%)	9 (0.3%)
24521	14 (0.6%)	0 (0.0%)	0 (0.0%)	14 (0.5%)
24800	71 (2.9%)	0 (0.0%)	0 (0.0%)	71 (2.3%)
25644	1 (0.0%)	0 (0.0%)	0 (0.0%)	1 (0.0%)
B	0 (0.0%)	0 (0.0%)	32 (33.0%)	32 (1.0%)
C	0 (0.0%)	0 (0.0%)	2 (2.1%)	2 (0.1%)
F	0 (0.0%)	0 (0.0%)	8 (8.2%)	8 (0.3%)
G	0 (0.0%)	0 (0.0%)	39 (40.2%)	39 (1.3%)
H	0 (0.0%)	0 (0.0%)	12 (12.4%)	12 (0.4%)
J	0 (0.0%)	0 (0.0%)	4 (4.1%)	4 (0.1%)

Pearson’s Chi-squared test

CAN <- bound3 %>% filter(COUNTRY == "CAN")
UK <- bound3 %>% filter(COUNTRY == "UK")
US <- bound3 %>% filter(COUNTRY == "US")

tab1 <- tableby(DTYPE ~ 
              factor(TX_YR),
                data=CAN, test=TRUE, total=TRUE, 
                numeric.stats=c("medianq1q3"), numeric.test="kwt", cat.test="chisq")
summary(tab1, title='Recipient characteristics', pfootnote=TRUE, digits = 2)

Recipient characteristics
	DBD (N=2770)	DCD (N=300)	LDLT (N=556)	Total (N=3626)	p value
factor(TX_YR)					< 0.001¹
2008	211 (7.6%)	16 (5.3%)	57 (10.3%)	284 (7.8%)
2009	216 (7.8%)	10 (3.3%)	47 (8.5%)	273 (7.5%)
2010	237 (8.6%)	12 (4.0%)	49 (8.8%)	298 (8.2%)
2011	252 (9.1%)	19 (6.3%)	47 (8.5%)	318 (8.8%)
2012	254 (9.2%)	20 (6.7%)	62 (11.2%)	336 (9.3%)
2013	234 (8.4%)	21 (7.0%)	49 (8.8%)	304 (8.4%)
2014	237 (8.6%)	44 (14.7%)	58 (10.4%)	339 (9.3%)
2015	249 (9.0%)	44 (14.7%)	54 (9.7%)	347 (9.6%)
2016	316 (11.4%)	29 (9.7%)	43 (7.7%)	388 (10.7%)
2017	298 (10.8%)	50 (16.7%)	43 (7.7%)	391 (10.8%)
2018	266 (9.6%)	35 (11.7%)	47 (8.5%)	348 (9.6%)

Pearson’s Chi-squared test

tab1 <- tableby(DTYPE ~ 
              factor(TX_YR),
                data=UK, test=TRUE, total=TRUE, 
                numeric.stats=c("medianq1q3"), numeric.test="kwt", cat.test="chisq")
summary(tab1, title='Recipient characteristics', pfootnote=TRUE, digits = 2)

Recipient characteristics
	DBD (N=4937)	DCD (N=1561)	LDLT (N=97)	Total (N=6595)	p value
factor(TX_YR)					< 0.001¹
2008	391 (7.9%)	74 (4.7%)	8 (8.2%)	473 (7.2%)
2009	374 (7.6%)	75 (4.8%)	10 (10.3%)	459 (7.0%)
2010	380 (7.7%)	100 (6.4%)	6 (6.2%)	486 (7.4%)
2011	394 (8.0%)	119 (7.6%)	10 (10.3%)	523 (7.9%)
2012	415 (8.4%)	138 (8.8%)	10 (10.3%)	563 (8.5%)
2013	466 (9.4%)	141 (9.0%)	11 (11.3%)	618 (9.4%)
2014	502 (10.2%)	162 (10.4%)	8 (8.2%)	672 (10.2%)
2015	422 (8.5%)	188 (12.0%)	15 (15.5%)	625 (9.5%)
2016	465 (9.4%)	202 (12.9%)	10 (10.3%)	677 (10.3%)
2017	544 (11.0%)	178 (11.4%)	6 (6.2%)	728 (11.0%)
2018	584 (11.8%)	184 (11.8%)	3 (3.1%)	771 (11.7%)

Pearson’s Chi-squared test

tab1 <- tableby(DTYPE ~ 
              factor(TX_YR),
                data=US, test=TRUE, total=TRUE, 
                numeric.stats=c("medianq1q3"), numeric.test="kwt", cat.test="chisq")
summary(tab1, title='Recipient characteristics', pfootnote=TRUE, digits = 2)

Recipient characteristics
	DBD (N=54151)	DCD (N=3601)	LDLT (N=2433)	Total (N=60185)	p value
factor(TX_YR)					< 0.001¹
2008	4501 (8.3%)	257 (7.1%)	162 (6.7%)	4920 (8.2%)
2009	4554 (8.4%)	256 (7.1%)	148 (6.1%)	4958 (8.2%)
2010	4516 (8.3%)	244 (6.8%)	196 (8.1%)	4956 (8.2%)
2011	4626 (8.5%)	244 (6.8%)	176 (7.2%)	5046 (8.4%)
2012	4555 (8.4%)	243 (6.7%)	181 (7.4%)	4979 (8.3%)
2013	4625 (8.5%)	287 (8.0%)	200 (8.2%)	5112 (8.5%)
2014	4786 (8.8%)	325 (9.0%)	213 (8.8%)	5324 (8.8%)
2015	4977 (9.2%)	368 (10.2%)	267 (11.0%)	5612 (9.3%)
2016	5522 (10.2%)	402 (11.2%)	274 (11.3%)	6198 (10.3%)
2017	5659 (10.5%)	478 (13.3%)	288 (11.8%)	6425 (10.7%)
2018	5830 (10.8%)	497 (13.8%)	328 (13.5%)	6655 (11.1%)

Pearson’s Chi-squared test

bound3 %>% select(DTYPE, DTYPE_DDLTandLDLT, GRAFT_TYPE) %>% View()

## Warning in system2("/usr/bin/otool", c("-L", shQuote(DSO)), stdout = TRUE):
## running command ''/usr/bin/otool' -L '/Library/Frameworks/R.framework/Resources/
## modules/R_de.so'' had status 1

Summaries for survival

#Overall
summary(tableby(COUNTRY ~ Surv(PSURV_years, PCENS), data=bound2, times=1:10, surv.stats=c("NeventsSurv", "NriskSurv")))

	US (N=2433)	CAN (N=556)	UK (N=97)	Total (N=3086)	p value
Surv(PSURV_years, PCENS)					0.001
time = 1	163 (93.0)	23 (95.8)	8 (91.4)	194 (93.5)
time = 2	221 (90.0)	32 (93.9)	9 (90.0)	262 (90.8)
time = 3	252 (88.1)	40 (92.1)	9 (90.0)	301 (88.9)
time = 4	284 (85.8)	44 (91.1)	10 (88.0)	338 (86.9)
time = 5	309 (83.5)	48 (89.9)	11 (85.4)	368 (84.8)
time = 6	333 (80.8)	52 (88.4)	11 (85.4)	396 (82.5)
time = 7	347 (78.9)	53 (87.9)	12 (80.0)	412 (80.8)
time = 8	351 (78.2)	59 (84.4)	12 (80.0)	422 (79.5)
time = 9	362 (75.3)	61 (82.7)	13 (66.7)	436 (76.8)
time = 10	373 (71.1)	61 (82.7)	13 (66.7)	447 (73.5)
time = 1	1979 (93.0)	484 (95.8)	79 (91.4)	2542 (93.5)
time = 2	1588 (90.0)	434 (93.9)	65 (90.0)	2087 (90.8)
time = 3	1313 (88.1)	384 (92.1)	50 (90.0)	1747 (88.9)
time = 4	1045 (85.8)	329 (91.1)	37 (88.0)	1411 (86.9)
time = 5	814 (83.5)	272 (89.9)	27 (85.4)	1113 (84.8)
time = 6	627 (80.8)	225 (88.4)	19 (85.4)	871 (82.5)
time = 7	485 (78.9)	168 (87.9)	12 (80.0)	665 (80.8)
time = 8	371 (78.2)	123 (84.4)	6 (80.0)	500 (79.5)
time = 9	229 (75.3)	82 (82.7)	4 (66.7)	315 (76.8)
time = 10	132 (71.1)	43 (82.7)	2 (66.7)	177 (73.5)

#1 year
summary(tableby(COUNTRY ~ Surv(PSURV_1year, PCENS_1year), data=bound2, times=1:1, surv.stats=c("NeventsSurv", "NriskSurv")))

	US (N=2433)	CAN (N=556)	UK (N=97)	Total (N=3086)	p value
Surv(PSURV_1year, PCENS_1year)					0.057
time = 1	163 (93.0)	23 (95.8)	8 (91.4)	194 (93.5)
time = 1	1979 (93.0)	484 (95.8)	79 (91.4)	2542 (93.5)

#3 year
summary(tableby(COUNTRY ~ Surv(PSURV_3year, PCENS_3year), data=bound2, times=1:3, surv.stats=c("NeventsSurv", "NriskSurv")))

	US (N=2433)	CAN (N=556)	UK (N=97)	Total (N=3086)	p value
Surv(PSURV_3year, PCENS_3year)					0.038
time = 1	163 (93.0)	23 (95.8)	8 (91.4)	194 (93.5)
time = 2	221 (90.0)	32 (93.9)	9 (90.0)	262 (90.8)
time = 3	252 (88.1)	40 (92.1)	9 (90.0)	301 (88.9)
time = 1	1979 (93.0)	484 (95.8)	79 (91.4)	2542 (93.5)
time = 2	1588 (90.0)	434 (93.9)	65 (90.0)	2087 (90.8)
time = 3	1313 (88.1)	384 (92.1)	50 (90.0)	1747 (88.9)

#5 year
summary(tableby(COUNTRY ~ Surv(PSURV_5year, PCENS_5year), data=bound2, times=1:5, surv.stats=c("NeventsSurv", "NriskSurv")))

	US (N=2433)	CAN (N=556)	UK (N=97)	Total (N=3086)	p value
Surv(PSURV_5year, PCENS_5year)					0.005
time = 1	163 (93.0)	23 (95.8)	8 (91.4)	194 (93.5)
time = 2	221 (90.0)	32 (93.9)	9 (90.0)	262 (90.8)
time = 3	252 (88.1)	40 (92.1)	9 (90.0)	301 (88.9)
time = 4	284 (85.8)	44 (91.1)	10 (88.0)	338 (86.9)
time = 5	309 (83.5)	48 (89.9)	11 (85.4)	368 (84.8)
time = 1	1979 (93.0)	484 (95.8)	79 (91.4)	2542 (93.5)
time = 2	1588 (90.0)	434 (93.9)	65 (90.0)	2087 (90.8)
time = 3	1313 (88.1)	384 (92.1)	50 (90.0)	1747 (88.9)
time = 4	1045 (85.8)	329 (91.1)	37 (88.0)	1411 (86.9)
time = 5	814 (83.5)	272 (89.9)	27 (85.4)	1113 (84.8)

#10 year
summary(tableby(COUNTRY ~ Surv(PSURV_10year, PCENS_10year), data=bound2, times=1:10, surv.stats=c("NeventsSurv", "NriskSurv")))

	US (N=2433)	CAN (N=556)	UK (N=97)	Total (N=3086)	p value
Surv(PSURV_10year, PCENS_10year)					0.001
time = 1	163 (93.0)	23 (95.8)	8 (91.4)	194 (93.5)
time = 2	221 (90.0)	32 (93.9)	9 (90.0)	262 (90.8)
time = 3	252 (88.1)	40 (92.1)	9 (90.0)	301 (88.9)
time = 4	284 (85.8)	44 (91.1)	10 (88.0)	338 (86.9)
time = 5	309 (83.5)	48 (89.9)	11 (85.4)	368 (84.8)
time = 6	333 (80.8)	52 (88.4)	11 (85.4)	396 (82.5)
time = 7	347 (78.9)	53 (87.9)	12 (80.0)	412 (80.8)
time = 8	351 (78.2)	59 (84.4)	12 (80.0)	422 (79.5)
time = 9	362 (75.3)	61 (82.7)	13 (66.7)	436 (76.8)
time = 10	373 (71.1)	61 (82.7)	13 (66.7)	447 (73.5)
time = 1	1979 (93.0)	484 (95.8)	79 (91.4)	2542 (93.5)
time = 2	1588 (90.0)	434 (93.9)	65 (90.0)	2087 (90.8)
time = 3	1313 (88.1)	384 (92.1)	50 (90.0)	1747 (88.9)
time = 4	1045 (85.8)	329 (91.1)	37 (88.0)	1411 (86.9)
time = 5	814 (83.5)	272 (89.9)	27 (85.4)	1113 (84.8)
time = 6	627 (80.8)	225 (88.4)	19 (85.4)	871 (82.5)
time = 7	485 (78.9)	168 (87.9)	12 (80.0)	665 (80.8)
time = 8	371 (78.2)	123 (84.4)	6 (80.0)	500 (79.5)
time = 9	229 (75.3)	82 (82.7)	4 (66.7)	315 (76.8)
time = 10	132 (71.1)	43 (82.7)	2 (66.7)	177 (73.5)

#Median
survfit(Surv(PSURV_years, PCENS) ~ COUNTRY, data = bound2)

Call: survfit(formula = Surv(PSURV_years, PCENS) ~ COUNTRY, data = bound2)

           n events median 0.95LCL 0.95UCL

COUNTRY=US 2433 377 NA NA NA COUNTRY=CAN 556 61 NA NA NA COUNTRY=UK 97 13 NA 8.13 NA

summary(survfit(Surv(PSURV_years, PCENS)~COUNTRY, data=bound2), times=1:10)

Call: survfit(formula = Surv(PSURV_years, PCENS) ~ COUNTRY, data = bound2)

            COUNTRY=US

time n.risk n.event survival std.err lower 95% CI upper 95% CI 1 1979 163 0.930 0.00528 0.920 0.941 2 1588 58 0.900 0.00640 0.888 0.913 3 1313 31 0.881 0.00712 0.868 0.895 4 1045 32 0.858 0.00805 0.842 0.874 5 814 25 0.835 0.00902 0.818 0.853 6 627 24 0.808 0.01027 0.788 0.829 7 485 14 0.789 0.01127 0.767 0.811 8 371 4 0.782 0.01173 0.759 0.805 9 229 11 0.753 0.01412 0.726 0.781 10 132 11 0.711 0.01823 0.676 0.748

            COUNTRY=CAN

time n.risk n.event survival std.err lower 95% CI upper 95% CI 1 484 23 0.958 0.00867 0.941 0.975 2 434 9 0.939 0.01045 0.919 0.960 3 384 8 0.921 0.01204 0.898 0.945 4 329 4 0.911 0.01300 0.886 0.936 5 272 4 0.899 0.01415 0.871 0.927 6 225 4 0.884 0.01571 0.854 0.915 7 168 1 0.879 0.01646 0.847 0.912 8 123 6 0.844 0.02110 0.804 0.886 9 82 2 0.827 0.02385 0.782 0.875 10 43 0 0.827 0.02385 0.782 0.875

            COUNTRY=UK

time n.risk n.event survival std.err lower 95% CI upper 95% CI 1 79 8 0.914 0.0293 0.858 0.973 2 65 1 0.900 0.0317 0.840 0.965 3 50 0 0.900 0.0317 0.840 0.965 4 37 1 0.880 0.0368 0.811 0.955 5 27 1 0.854 0.0443 0.771 0.945 6 19 0 0.854 0.0443 0.771 0.945 7 12 1 0.800 0.0663 0.680 0.941 8 6 0 0.800 0.0663 0.680 0.941 9 4 1 0.667 0.1337 0.450 0.988 10 2 0 0.667 0.1337 0.450 0.988

fit4 <- pairwise_survdiff(Surv(PSURV_years, PCENS) ~ COUNTRY , data = bound2)
fit4

Pairwise comparisons using Log-Rank test

data: bound2 and COUNTRY

US      CAN

CAN 0.00066 -
UK 0.86368 0.19571

P value adjustment method: BH

symnum(fit4$p.value, cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 0.1, 1),
   symbols = c("****", "***", "**", "*", "+", " "),
   abbr.colnames = FALSE, na = "")

US  CAN

CAN ***
UK
attr(,“legend”) [1] 0 ‘****’ 0.0001 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘+’ 0.1 ’ ’ 1 \t ## NA: ’’

Summaries for graft survival

#Overall
summary(tableby(COUNTRY ~ Surv(GSURV_years, GCENS), data=bound2, times=1:10, surv.stats=c("NeventsSurv", "NriskSurv")))

	US (N=2433)	CAN (N=556)	UK (N=97)	Total (N=3086)	p value
Surv(GSURV_years, GCENS)					< 0.001
time = 1	282 (88.4)	39 (92.8)	11 (88.2)	332 (89.1)
time = 2	353 (84.9)	49 (90.8)	12 (86.7)	414 (86.0)
time = 3	387 (83.0)	55 (89.4)	12 (86.7)	454 (84.2)
time = 4	425 (80.3)	61 (87.8)	12 (86.7)	498 (81.8)
time = 5	455 (77.8)	62 (87.4)	12 (86.7)	529 (79.8)
time = 6	481 (75.0)	64 (86.7)	12 (86.7)	557 (77.5)
time = 7	497 (73.0)	66 (85.8)	13 (81.3)	576 (75.7)
time = 8	505 (71.6)	69 (83.8)	13 (81.3)	587 (74.2)
time = 9	517 (68.8)	72 (81.1)	13 (81.3)	602 (71.4)
time = 10	529 (64.6)	72 (81.1)	13 (81.3)	614 (68.0)
time = 1	1979 (88.4)	465 (92.8)	71 (88.2)	2515 (89.1)
time = 2	1587 (84.9)	412 (90.8)	56 (86.7)	2055 (86.0)
time = 3	1313 (83.0)	362 (89.4)	42 (86.7)	1717 (84.2)
time = 4	1044 (80.3)	305 (87.8)	32 (86.7)	1381 (81.8)
time = 5	815 (77.8)	250 (87.4)	23 (86.7)	1088 (79.8)
time = 6	627 (75.0)	210 (86.7)	17 (86.7)	854 (77.5)
time = 7	485 (73.0)	156 (85.8)	9 (81.3)	650 (75.7)
time = 8	371 (71.6)	114 (83.8)	6 (81.3)	491 (74.2)
time = 9	228 (68.8)	74 (81.1)	4 (81.3)	306 (71.4)
time = 10	131 (64.6)	38 (81.1)	2 (81.3)	171 (68.0)

#1 year
summary(tableby(COUNTRY ~ Surv(GSURV_1year, GCENS_1year), data=bound2, times=1:1, surv.stats=c("NeventsSurv", "NriskSurv")))

	US (N=2433)	CAN (N=556)	UK (N=97)	Total (N=3086)	p value
Surv(GSURV_1year, GCENS_1year)					0.015
time = 1	282 (88.4)	39 (92.8)	11 (88.2)	332 (89.1)
time = 1	1979 (88.4)	465 (92.8)	71 (88.2)	2515 (89.1)

#3 year
summary(tableby(COUNTRY ~ Surv(GSURV_3year, GCENS_3year), data=bound2, times=1:3, surv.stats=c("NeventsSurv", "NriskSurv")))

	US (N=2433)	CAN (N=556)	UK (N=97)	Total (N=3086)	p value
Surv(GSURV_3year, GCENS_3year)					0.002
time = 1	282 (88.4)	39 (92.8)	11 (88.2)	332 (89.1)
time = 2	353 (84.9)	49 (90.8)	12 (86.7)	414 (86.0)
time = 3	387 (83.0)	55 (89.4)	12 (86.7)	454 (84.2)
time = 1	1979 (88.4)	465 (92.8)	71 (88.2)	2515 (89.1)
time = 2	1587 (84.9)	412 (90.8)	56 (86.7)	2055 (86.0)
time = 3	1313 (83.0)	362 (89.4)	42 (86.7)	1717 (84.2)

#5 year
summary(tableby(COUNTRY ~ Surv(GSURV_5year, GCENS_5year), data=bound2, times=1:5, surv.stats=c("NeventsSurv", "NriskSurv")))

	US (N=2433)	CAN (N=556)	UK (N=97)	Total (N=3086)	p value
Surv(GSURV_5year, GCENS_5year)					< 0.001
time = 1	282 (88.4)	39 (92.8)	11 (88.2)	332 (89.1)
time = 2	353 (84.9)	49 (90.8)	12 (86.7)	414 (86.0)
time = 3	387 (83.0)	55 (89.4)	12 (86.7)	454 (84.2)
time = 4	425 (80.3)	61 (87.8)	12 (86.7)	498 (81.8)
time = 5	455 (77.8)	62 (87.4)	12 (86.7)	529 (79.8)
time = 1	1979 (88.4)	465 (92.8)	71 (88.2)	2515 (89.1)
time = 2	1587 (84.9)	412 (90.8)	56 (86.7)	2055 (86.0)
time = 3	1313 (83.0)	362 (89.4)	42 (86.7)	1717 (84.2)
time = 4	1044 (80.3)	305 (87.8)	32 (86.7)	1381 (81.8)
time = 5	815 (77.8)	250 (87.4)	23 (86.7)	1088 (79.8)

#10 year
summary(tableby(COUNTRY ~ Surv(GSURV_10year, GCENS_10year), data=bound2, times=1:10, surv.stats=c("NeventsSurv", "NriskSurv")))

	US (N=2433)	CAN (N=556)	UK (N=97)	Total (N=3086)	p value
Surv(GSURV_10year, GCENS_10year)					< 0.001
time = 1	282 (88.4)	39 (92.8)	11 (88.2)	332 (89.1)
time = 2	353 (84.9)	49 (90.8)	12 (86.7)	414 (86.0)
time = 3	387 (83.0)	55 (89.4)	12 (86.7)	454 (84.2)
time = 4	425 (80.3)	61 (87.8)	12 (86.7)	498 (81.8)
time = 5	455 (77.8)	62 (87.4)	12 (86.7)	529 (79.8)
time = 6	481 (75.0)	64 (86.7)	12 (86.7)	557 (77.5)
time = 7	497 (73.0)	66 (85.8)	13 (81.3)	576 (75.7)
time = 8	505 (71.6)	69 (83.8)	13 (81.3)	587 (74.2)
time = 9	517 (68.8)	72 (81.1)	13 (81.3)	602 (71.4)
time = 10	529 (64.6)	72 (81.1)	13 (81.3)	614 (68.0)
time = 1	1979 (88.4)	465 (92.8)	71 (88.2)	2515 (89.1)
time = 2	1587 (84.9)	412 (90.8)	56 (86.7)	2055 (86.0)
time = 3	1313 (83.0)	362 (89.4)	42 (86.7)	1717 (84.2)
time = 4	1044 (80.3)	305 (87.8)	32 (86.7)	1381 (81.8)
time = 5	815 (77.8)	250 (87.4)	23 (86.7)	1088 (79.8)
time = 6	627 (75.0)	210 (86.7)	17 (86.7)	854 (77.5)
time = 7	485 (73.0)	156 (85.8)	9 (81.3)	650 (75.7)
time = 8	371 (71.6)	114 (83.8)	6 (81.3)	491 (74.2)
time = 9	228 (68.8)	74 (81.1)	4 (81.3)	306 (71.4)
time = 10	131 (64.6)	38 (81.1)	2 (81.3)	171 (68.0)

#Median
survfit(Surv(GSURV_years, GCENS) ~ COUNTRY, data = bound2)

Call: survfit(formula = Surv(GSURV_years, GCENS) ~ COUNTRY, data = bound2)

           n events median 0.95LCL 0.95UCL

COUNTRY=US 2433 537 NA 11.1 NA COUNTRY=CAN 556 72 NA NA NA COUNTRY=UK 97 13 NA NA NA

summary(survfit(Surv(GSURV_years, GCENS)~COUNTRY, data=bound2), times=1:10)

Call: survfit(formula = Surv(GSURV_years, GCENS) ~ COUNTRY, data = bound2)

            COUNTRY=US

time n.risk n.event survival std.err lower 95% CI upper 95% CI 1 1979 282 0.884 0.00652 0.871 0.896 2 1587 71 0.849 0.00744 0.835 0.864 3 1313 34 0.830 0.00799 0.814 0.845 4 1044 38 0.803 0.00881 0.786 0.821 5 815 30 0.778 0.00969 0.759 0.797 6 627 26 0.750 0.01072 0.730 0.772 7 485 16 0.730 0.01161 0.707 0.753 8 371 8 0.716 0.01235 0.692 0.741 9 228 12 0.688 0.01430 0.661 0.717 10 131 12 0.646 0.01795 0.611 0.682

            COUNTRY=CAN

time n.risk n.event survival std.err lower 95% CI upper 95% CI 1 465 39 0.928 0.0111 0.907 0.950 2 412 10 0.908 0.0126 0.883 0.933 3 362 6 0.894 0.0136 0.868 0.921 4 305 6 0.878 0.0149 0.849 0.907 5 250 1 0.874 0.0152 0.845 0.905 6 210 2 0.867 0.0160 0.836 0.899 7 156 2 0.858 0.0171 0.825 0.892 8 114 3 0.838 0.0200 0.800 0.879 9 74 3 0.811 0.0249 0.764 0.861 10 38 0 0.811 0.0249 0.764 0.861

            COUNTRY=UK

time n.risk n.event survival std.err lower 95% CI upper 95% CI 1 71 11 0.882 0.0336 0.818 0.950 2 56 1 0.867 0.0360 0.800 0.941 3 42 0 0.867 0.0360 0.800 0.941 4 32 0 0.867 0.0360 0.800 0.941 5 23 0 0.867 0.0360 0.800 0.941 6 17 0 0.867 0.0360 0.800 0.941 7 9 1 0.813 0.0624 0.700 0.945 8 6 0 0.813 0.0624 0.700 0.945 9 4 0 0.813 0.0624 0.700 0.945 10 2 0 0.813 0.0624 0.700 0.945

fit4 <- pairwise_survdiff(Surv(GSURV_years, GCENS) ~ COUNTRY , data = bound2)
fit4

Pairwise comparisons using Log-Rank test

data: bound2 and COUNTRY

US        CAN

CAN 0.0000018 -
UK 0.35 0.45

P value adjustment method: BH

symnum(fit4$p.value, cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 0.1, 1),
   symbols = c("****", "***", "**", "*", "+", " "),
   abbr.colnames = FALSE, na = "")

US   CAN

CAN ****
UK
attr(,“legend”) [1] 0 ‘****’ 0.0001 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘+’ 0.1 ’ ’ 1 \t ## NA: ’’

OS overall Figure 3

levels(bound2$COUNTRY)

[1] “US” “CAN” “UK”

fit1 <- survfit(Surv(PSURV_years, PCENS) ~ COUNTRY, data = bound2)

ggsurv1 <- ggsurvplot(
          fit1,               
          data = bound2,          
          risk.table = TRUE,     
          pval = TRUE,
          pval.size = 6,
          pval.coord = c(0, 0.8),
          conf.int = F,         
          xlim = c(0,5),
          ylim = c(0.7, 1.00),
          xlab = "Years from LT",
          palette = "jama",
          break.time.by = 1,  
          title = "",
          ggtheme = theme_test(base_size=28, base_family = "Helvetica"),
          fontsize=8,
          risk.table.height = 0.25,
          tables.x.text = FALSE,
          legend.title= "",
          legend.labs =
           c("US", "Canada", "UK"),
        )

## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.

ggsurv1$table <- ggsurv1$table + labs(x = NULL, y = NULL) + theme(plot.title = element_text(size = 28))

ggsurv1$plot <- ggsurv1$plot+ 
              ggplot2::annotate("text", 
                                x = 2.5, y = 0.7, # x and y coordinates of the text
                                label = "Adjusted HR (ref: US, Canada HR 0.53, 95% CI 0.25-1.14, UK HR 1.34, 95% CI 0.71-2.53)", size = 6) 
ggsurv1

## Warning: Removed 72 row(s) containing missing values (geom_path).

## Warning: Removed 69 rows containing missing values (geom_point).

## Warning: Removed 72 row(s) containing missing values (geom_path).

## Warning: Removed 69 rows containing missing values (geom_point).

fit4 <- pairwise_survdiff(Surv(PSURV_years, PCENS) ~ COUNTRY, data = bound2)
fit4

Pairwise comparisons using Log-Rank test

data: bound2 and COUNTRY

US      CAN

CAN 0.00066 -
UK 0.86368 0.19571

P value adjustment method: BH

symnum(fit4$p.value, cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 0.1, 1),
   symbols = c("****", "***", "**", "*", "+", " "),
   abbr.colnames = FALSE, na = "")

US  CAN

CAN ***
UK
attr(,“legend”) [1] 0 ‘****’ 0.0001 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘+’ 0.1 ’ ’ 1 \t ## NA: ’’

Graft survival overall Figure S1

levels(bound2$COUNTRY)

[1] “US” “CAN” “UK”

fit1 <- survfit(Surv(GSURV_years, GCENS) ~ COUNTRY, data = bound2)

ggsurv1 <- ggsurvplot(
          fit1,               
          data = bound2,          
          risk.table = TRUE,     
          pval = TRUE,
          pval.size = 6,
          pval.coord = c(0, 0.8),
          conf.int = F,         
          xlim = c(0,5),
          ylim = c(0.7, 1.00),
          xlab = "Years from LT",
          palette = "jama",
          break.time.by = 1,  
          title = "",
          ggtheme = theme_test(base_size=28, base_family = "Helvetica"),
          fontsize=8,
          risk.table.height = 0.25,
          tables.x.text = FALSE,
          legend.title= "",
          legend.labs =
           c("US", "Canada", "UK"),
        )

## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.

ggsurv1$table <- ggsurv1$table + labs(x = NULL, y = NULL) + theme(plot.title = element_text(size = 28))

ggsurv1$plot <- ggsurv1$plot+ 
              ggplot2::annotate("text", 
                                x = 2.5, y = 0.7, # x and y coordinates of the text
                                label = "Adjusted HR (ref: US, Canada HR 0.54, 95% CI 0.29-1.01, UK HR 0.91, 95% CI 0.49-1.70)", size = 6) 
ggsurv1

## Warning: Removed 219 row(s) containing missing values (geom_path).

## Warning: Removed 199 rows containing missing values (geom_point).

## Warning: Removed 219 row(s) containing missing values (geom_path).

## Warning: Removed 199 rows containing missing values (geom_point).

fit4 <- pairwise_survdiff(Surv(GSURV_years, GCENS) ~ COUNTRY, data = bound2)
fit4

Pairwise comparisons using Log-Rank test

data: bound2 and COUNTRY

US        CAN

CAN 0.0000018 -
UK 0.35 0.45

P value adjustment method: BH

symnum(fit4$p.value, cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 0.1, 1),
   symbols = c("****", "***", "**", "*", "+", " "),
   abbr.colnames = FALSE, na = "")

US   CAN

CAN ****
UK
attr(,“legend”) [1] 0 ‘****’ 0.0001 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘+’ 0.1 ’ ’ 1 \t ## NA: ’’

OS overall by percentile

levels(bound2$COUNTRY)

[1] “US” “CAN” “UK”

fit1 <- survfit(Surv(PSURV_years, PCENS) ~ Centerpercentile, data = bound2)

ggsurv1 <- ggsurvplot(
          fit1,               
          data = bound2,          
          risk.table = TRUE,     
          pval = TRUE,
          pval.size = 6,
          pval.coord = c(0, 0.8),
          conf.int = F,         
          xlim = c(0,5),
          ylim = c(0.7, 1.00),
          xlab = "Years from LT",
          palette = "jama",
          break.time.by = 1,  
          title = "",
          ggtheme = theme_test(base_size=28, base_family = "Helvetica"),
          fontsize=8,
          risk.table.height = 0.25,
          tables.x.text = FALSE,
          legend.title= "",
          legend.labs =
           c("<25%ile", "25%-75%ile", ">75%ile"),
        )

## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.

ggsurv1

## Warning: Removed 28 row(s) containing missing values (geom_path).

## Warning: Removed 27 rows containing missing values (geom_point).

## Warning: Removed 28 row(s) containing missing values (geom_path).

## Warning: Removed 27 rows containing missing values (geom_point).

OS overall DDLT vs LDLT US Figure 4a

US <- bound3 %>% filter(COUNTRY == "US")
fit1 <- survfit(Surv(PSURV_years, PCENS) ~ LDLT, data = US)

ggsurv1 <- ggsurvplot(
           fit1,               
           data = US,          
           risk.table = TRUE,     
           pval = TRUE,
           pval.coord = c(0, 0.8),
           pval.size = 6,
           conf.int = F,         
           xlim = c(0,5),
           ylim = c(0.5, 1.00),
           xlab = "Years from LT",
           palette = "jama",
           break.time.by = 1,  
           title = "United States",
      # subtitle = "with 95% confidence intervals",
               ggtheme = theme_test(base_size=28, base_family = "Helvetica"),
         # risk.table.y.text.col = F,
         fontsize=8,
         risk.table.height = 0.25,
        #  risk.table.y.text = T,
        #  surv.median.line = "hv", 
           legend.title= "",
      tables.x.text = FALSE,
          legend.labs =
           c("DDLT", "LDLT")
        )

## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.

ggsurv1

ggsurv1$plot <- ggsurv1$plot+ 
              ggplot2::annotate("text", 
                                x = 2.5, y = 0.5, # x and y coordinates of the text
                                label = "Adjusted HR (ref: DDLT, HR LDLT 1.03, 95% CI 0.92-1.16)", size = 6)
ggsurv1

#Overall
summary(tableby(LDLT ~ Surv(PSURV_years, PCENS), data=US, times=1:10, surv.stats=c("NeventsSurv", "NriskSurv")))

	0 (N=57752)	1 (N=2433)	Total (N=60185)	p value
Surv(PSURV_years, PCENS)				< 0.001
time = 1	4795 (91.5)	163 (93.0)	4958 (91.6)
time = 2	6758 (87.6)	221 (90.0)	6979 (87.7)
time = 3	8194 (84.2)	252 (88.1)	8446 (84.4)
time = 4	9275 (81.2)	284 (85.8)	9559 (81.4)
time = 5	10208 (78.2)	309 (83.5)	10517 (78.4)
time = 6	10935 (75.3)	333 (80.8)	11268 (75.5)
time = 7	11509 (72.5)	347 (78.9)	11856 (72.8)
time = 8	11947 (69.7)	351 (78.2)	12298 (70.0)
time = 9	12284 (66.8)	362 (75.3)	12646 (67.1)
time = 10	12543 (63.3)	373 (71.1)	12916 (63.5)
time = 1	48435 (91.5)	1979 (93.0)	50414 (91.6)
time = 2	40076 (87.6)	1588 (90.0)	41664 (87.7)
time = 3	32852 (84.2)	1313 (88.1)	34165 (84.4)
time = 4	26844 (81.2)	1045 (85.8)	27889 (81.4)
time = 5	21520 (78.2)	814 (83.5)	22334 (78.4)
time = 6	16843 (75.3)	627 (80.8)	17470 (75.5)
time = 7	12795 (72.5)	485 (78.9)	13280 (72.8)
time = 8	9153 (69.7)	371 (78.2)	9524 (70.0)
time = 9	5899 (66.8)	229 (75.3)	6128 (67.1)
time = 10	3126 (63.3)	132 (71.1)	3258 (63.5)

#1 year
summary(tableby(LDLT ~ Surv(PSURV_1year, PCENS_1year), data=US, times=1:1, surv.stats=c("NeventsSurv", "NriskSurv")))

	0 (N=57752)	1 (N=2433)	Total (N=60185)	p value
Surv(PSURV_1year, PCENS_1year)				0.010
time = 1	4795 (91.5)	163 (93.0)	4958 (91.6)
time = 1	48435 (91.5)	1979 (93.0)	50414 (91.6)

#3 year
summary(tableby(LDLT ~ Surv(PSURV_3year, PCENS_3year), data=US, times=1:3, surv.stats=c("NeventsSurv", "NriskSurv")))

	0 (N=57752)	1 (N=2433)	Total (N=60185)	p value
Surv(PSURV_3year, PCENS_3year)				< 0.001
time = 1	4795 (91.5)	163 (93.0)	4958 (91.6)
time = 2	6758 (87.6)	221 (90.0)	6979 (87.7)
time = 3	8194 (84.2)	252 (88.1)	8446 (84.4)
time = 1	48435 (91.5)	1979 (93.0)	50414 (91.6)
time = 2	40076 (87.6)	1588 (90.0)	41664 (87.7)
time = 3	32852 (84.2)	1313 (88.1)	34165 (84.4)

#5 year
summary(tableby(LDLT ~ Surv(PSURV_5year, PCENS_5year), data=US, times=1:5, surv.stats=c("NeventsSurv", "NriskSurv")))

	0 (N=57752)	1 (N=2433)	Total (N=60185)	p value
Surv(PSURV_5year, PCENS_5year)				< 0.001
time = 1	4795 (91.5)	163 (93.0)	4958 (91.6)
time = 2	6758 (87.6)	221 (90.0)	6979 (87.7)
time = 3	8194 (84.2)	252 (88.1)	8446 (84.4)
time = 4	9275 (81.2)	284 (85.8)	9559 (81.4)
time = 5	10208 (78.2)	309 (83.5)	10517 (78.4)
time = 1	48435 (91.5)	1979 (93.0)	50414 (91.6)
time = 2	40076 (87.6)	1588 (90.0)	41664 (87.7)
time = 3	32852 (84.2)	1313 (88.1)	34165 (84.4)
time = 4	26844 (81.2)	1045 (85.8)	27889 (81.4)
time = 5	21520 (78.2)	814 (83.5)	22334 (78.4)

#10 year
summary(tableby(LDLT ~ Surv(PSURV_10year, PCENS_10year), data=US, times=1:10, surv.stats=c("NeventsSurv", "NriskSurv")))

	0 (N=57752)	1 (N=2433)	Total (N=60185)	p value
Surv(PSURV_10year, PCENS_10year)				< 0.001
time = 1	4795 (91.5)	163 (93.0)	4958 (91.6)
time = 2	6758 (87.6)	221 (90.0)	6979 (87.7)
time = 3	8194 (84.2)	252 (88.1)	8446 (84.4)
time = 4	9275 (81.2)	284 (85.8)	9559 (81.4)
time = 5	10208 (78.2)	309 (83.5)	10517 (78.4)
time = 6	10935 (75.3)	333 (80.8)	11268 (75.5)
time = 7	11509 (72.5)	347 (78.9)	11856 (72.8)
time = 8	11947 (69.7)	351 (78.2)	12298 (70.0)
time = 9	12284 (66.8)	362 (75.3)	12646 (67.1)
time = 10	12543 (63.3)	373 (71.1)	12916 (63.5)
time = 1	48435 (91.5)	1979 (93.0)	50414 (91.6)
time = 2	40076 (87.6)	1588 (90.0)	41664 (87.7)
time = 3	32852 (84.2)	1313 (88.1)	34165 (84.4)
time = 4	26844 (81.2)	1045 (85.8)	27889 (81.4)
time = 5	21520 (78.2)	814 (83.5)	22334 (78.4)
time = 6	16843 (75.3)	627 (80.8)	17470 (75.5)
time = 7	12795 (72.5)	485 (78.9)	13280 (72.8)
time = 8	9153 (69.7)	371 (78.2)	9524 (70.0)
time = 9	5899 (66.8)	229 (75.3)	6128 (67.1)
time = 10	3126 (63.3)	132 (71.1)	3258 (63.5)

#Median
survfit(Surv(PSURV_years, PCENS) ~ LDLT, data = US)

Call: survfit(formula = Surv(PSURV_years, PCENS) ~ LDLT, data = US)

       n events median 0.95LCL 0.95UCL

LDLT=0 57752 12708 NA 12 NA LDLT=1 2433 377 NA NA NA

summary(survfit(Surv(PSURV_years, PCENS)~LDLT, data=US), times=1:10)

Call: survfit(formula = Surv(PSURV_years, PCENS) ~ LDLT, data = US)

            LDLT=0

time n.risk n.event survival std.err lower 95% CI upper 95% CI 1 48435 4795 0.915 0.00117 0.913 0.918 2 40076 1963 0.876 0.00142 0.873 0.879 3 32852 1436 0.842 0.00162 0.839 0.845 4 26844 1081 0.812 0.00180 0.809 0.816 5 21520 933 0.782 0.00199 0.778 0.786 6 16843 727 0.753 0.00218 0.749 0.758 7 12795 574 0.725 0.00240 0.721 0.730 8 9153 438 0.697 0.00265 0.692 0.703 9 5899 337 0.668 0.00299 0.662 0.674 10 3126 259 0.633 0.00355 0.626 0.640

            LDLT=1

fit4 <- pairwise_survdiff(Surv(PSURV_years, PCENS) ~ LDLT , data = US)
fit4

Pairwise comparisons using Log-Rank test

data: US and LDLT

0
1 0.0000000041

P value adjustment method: BH

symnum(fit4$p.value, cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 0.1, 1),
   symbols = c("****", "***", "**", "*", "+", " "),
   abbr.colnames = FALSE, na = "")

0
1 **** attr(,“legend”) [1] 0 ‘****’ 0.0001 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘+’ 0.1 ’ ’ 1 OS overall DDLT vs LDLT CAN Fig 4b

CAN <- bound3 %>% filter(COUNTRY == "CAN")
fit1 <- survfit(Surv(PSURV_years, PCENS) ~ LDLT, data = CAN)

ggsurv1 <- ggsurvplot(
           fit1,               
           data = CAN,          
           risk.table = TRUE,     
           pval = TRUE,
           pval.coord = c(0, 0.8),
           pval.size = 6,
           conf.int = F,         
           xlim = c(0,5),
           ylim = c(0.5, 1.00),
           xlab = "Years from LT",
           palette = "jama",
           break.time.by = 1,  
           title = "Canada",
      # subtitle = "with 95% confidence intervals",
               ggtheme = theme_test(base_size=28, base_family = "Helvetica"),
         # risk.table.y.text.col = F,
         fontsize=8,
         risk.table.height = 0.25,
        #  risk.table.y.text = T,
        #  surv.median.line = "hv", 
           legend.title= "",
      tables.x.text = FALSE,
          legend.labs =
           c("DDLT", "LDLT")
        )

## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.

ggsurv1

ggsurv1$plot <- ggsurv1$plot+ 
              ggplot2::annotate("text", 
                                x = 2.5, y = 0.5, # x and y coordinates of the text
                                label = "Adjusted HR (ref: DDLT, HR LDLT 0.65, 95% CI 0.30-1.41)", size = 6)
ggsurv1

#Overall
summary(tableby(LDLT ~ Surv(PSURV_years, PCENS), data=CAN, times=1:10, surv.stats=c("NeventsSurv", "NriskSurv")))

	0 (N=3070)	1 (N=556)	Total (N=3626)	p value
Surv(PSURV_years, PCENS)				< 0.001
time = 1	250 (91.6)	23 (95.8)	273 (92.2)
time = 2	332 (88.4)	32 (93.9)	364 (89.3)
time = 3	378 (86.4)	40 (92.1)	418 (87.3)
time = 4	414 (84.5)	44 (91.1)	458 (85.5)
time = 5	444 (82.6)	48 (89.9)	492 (83.8)
time = 6	462 (81.3)	52 (88.4)	514 (82.4)
time = 7	486 (79.2)	53 (87.9)	539 (80.6)
time = 8	498 (77.6)	59 (84.4)	557 (78.7)
time = 9	504 (76.5)	61 (82.7)	565 (77.5)
time = 10	508 (75.4)	61 (82.7)	569 (76.5)
time = 1	2528 (91.6)	484 (95.8)	3012 (92.2)
time = 2	2124 (88.4)	434 (93.9)	2558 (89.3)
time = 3	1770 (86.4)	384 (92.1)	2154 (87.3)
time = 4	1495 (84.5)	329 (91.1)	1824 (85.5)
time = 5	1234 (82.6)	272 (89.9)	1506 (83.8)
time = 6	1018 (81.3)	225 (88.4)	1243 (82.4)
time = 7	763 (79.2)	168 (87.9)	931 (80.6)
time = 8	536 (77.6)	123 (84.4)	659 (78.7)
time = 9	350 (76.5)	82 (82.7)	432 (77.5)
time = 10	178 (75.4)	43 (82.7)	221 (76.5)

#1 year
summary(tableby(LDLT ~ Surv(PSURV_1year, PCENS_1year), data=CAN, times=1:1, surv.stats=c("NeventsSurv", "NriskSurv")))

	0 (N=3070)	1 (N=556)	Total (N=3626)	p value
Surv(PSURV_1year, PCENS_1year)				0.001
time = 1	250 (91.6)	23 (95.8)	273 (92.2)
time = 1	2528 (91.6)	484 (95.8)	3012 (92.2)

#3 year
summary(tableby(LDLT ~ Surv(PSURV_3year, PCENS_3year), data=CAN, times=1:3, surv.stats=c("NeventsSurv", "NriskSurv")))

	0 (N=3070)	1 (N=556)	Total (N=3626)	p value
Surv(PSURV_3year, PCENS_3year)				< 0.001
time = 1	250 (91.6)	23 (95.8)	273 (92.2)
time = 2	332 (88.4)	32 (93.9)	364 (89.3)
time = 3	378 (86.4)	40 (92.1)	418 (87.3)
time = 1	2528 (91.6)	484 (95.8)	3012 (92.2)
time = 2	2124 (88.4)	434 (93.9)	2558 (89.3)
time = 3	1770 (86.4)	384 (92.1)	2154 (87.3)

#5 year
summary(tableby(LDLT ~ Surv(PSURV_5year, PCENS_5year), data=CAN, times=1:5, surv.stats=c("NeventsSurv", "NriskSurv")))

	0 (N=3070)	1 (N=556)	Total (N=3626)	p value
Surv(PSURV_5year, PCENS_5year)				< 0.001
time = 1	250 (91.6)	23 (95.8)	273 (92.2)
time = 2	332 (88.4)	32 (93.9)	364 (89.3)
time = 3	378 (86.4)	40 (92.1)	418 (87.3)
time = 4	414 (84.5)	44 (91.1)	458 (85.5)
time = 5	444 (82.6)	48 (89.9)	492 (83.8)
time = 1	2528 (91.6)	484 (95.8)	3012 (92.2)
time = 2	2124 (88.4)	434 (93.9)	2558 (89.3)
time = 3	1770 (86.4)	384 (92.1)	2154 (87.3)
time = 4	1495 (84.5)	329 (91.1)	1824 (85.5)
time = 5	1234 (82.6)	272 (89.9)	1506 (83.8)

#10 year
summary(tableby(LDLT ~ Surv(PSURV_10year, PCENS_10year), data=CAN, times=1:10, surv.stats=c("NeventsSurv", "NriskSurv")))

	0 (N=3070)	1 (N=556)	Total (N=3626)	p value
Surv(PSURV_10year, PCENS_10year)				< 0.001
time = 1	250 (91.6)	23 (95.8)	273 (92.2)
time = 2	332 (88.4)	32 (93.9)	364 (89.3)
time = 3	378 (86.4)	40 (92.1)	418 (87.3)
time = 4	414 (84.5)	44 (91.1)	458 (85.5)
time = 5	444 (82.6)	48 (89.9)	492 (83.8)
time = 6	462 (81.3)	52 (88.4)	514 (82.4)
time = 7	486 (79.2)	53 (87.9)	539 (80.6)
time = 8	498 (77.6)	59 (84.4)	557 (78.7)
time = 9	504 (76.5)	61 (82.7)	565 (77.5)
time = 10	508 (75.4)	61 (82.7)	569 (76.5)
time = 1	2528 (91.6)	484 (95.8)	3012 (92.2)
time = 2	2124 (88.4)	434 (93.9)	2558 (89.3)
time = 3	1770 (86.4)	384 (92.1)	2154 (87.3)
time = 4	1495 (84.5)	329 (91.1)	1824 (85.5)
time = 5	1234 (82.6)	272 (89.9)	1506 (83.8)
time = 6	1018 (81.3)	225 (88.4)	1243 (82.4)
time = 7	763 (79.2)	168 (87.9)	931 (80.6)
time = 8	536 (77.6)	123 (84.4)	659 (78.7)
time = 9	350 (76.5)	82 (82.7)	432 (77.5)
time = 10	178 (75.4)	43 (82.7)	221 (76.5)

#Median
survfit(Surv(PSURV_years, PCENS) ~ LDLT, data = CAN)

Call: survfit(formula = Surv(PSURV_years, PCENS) ~ LDLT, data = CAN)

      n events median 0.95LCL 0.95UCL

LDLT=0 3070 509 NA NA NA LDLT=1 556 61 NA NA NA

summary(survfit(Surv(PSURV_years, PCENS)~LDLT, data=CAN), times=1:10)

Call: survfit(formula = Surv(PSURV_years, PCENS) ~ LDLT, data = CAN)

            LDLT=0

time n.risk n.event survival std.err lower 95% CI upper 95% CI 1 2528 250 0.916 0.00509 0.906 0.926 2 2124 82 0.884 0.00599 0.873 0.896 3 1770 46 0.864 0.00658 0.851 0.877 4 1495 36 0.845 0.00716 0.831 0.859 5 1234 30 0.826 0.00776 0.811 0.842 6 1018 18 0.813 0.00824 0.797 0.829 7 763 24 0.792 0.00911 0.774 0.810 8 536 12 0.776 0.00996 0.757 0.796 9 350 6 0.765 0.01079 0.745 0.787 10 178 4 0.754 0.01220 0.730 0.778

            LDLT=1

fit4 <- pairwise_survdiff(Surv(PSURV_years, PCENS) ~ LDLT , data = CAN)
fit4

Pairwise comparisons using Log-Rank test

data: CAN and LDLT

0
1 0.00017

P value adjustment method: BH

symnum(fit4$p.value, cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 0.1, 1),
   symbols = c("****", "***", "**", "*", "+", " "),
   abbr.colnames = FALSE, na = "")

0
1 *** attr(,“legend”) [1] 0 ‘****’ 0.0001 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘+’ 0.1 ’ ’ 1

OS overall DDLT vs LDLT UK Figure 4c

UK <- bound3 %>% filter(COUNTRY == "UK")
fit1 <- survfit(Surv(PSURV_years, PCENS) ~ LDLT, data = UK)

ggsurv1 <- ggsurvplot(
           fit1,               
           data = UK,          
           risk.table = TRUE,     
           pval = TRUE,
           pval.coord = c(0, 0.8),
           pval.size = 6,
           conf.int = F,         
           xlim = c(0,5),
           ylim = c(0.5, 1.00),
           xlab = "Years from LT",
           palette = "jama",
           break.time.by = 1,  
           title = "United Kingdom",
      # subtitle = "with 95% confidence intervals",
               ggtheme = theme_test(base_size=28, base_family = "Helvetica"),
         # risk.table.y.text.col = F,
         fontsize=8,
         risk.table.height = 0.25,
        #  risk.table.y.text = T,
        #  surv.median.line = "hv", 
           legend.title= "",
      tables.x.text = FALSE,
          legend.labs =
           c("DDLT", "LDLT")
        )

## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.

ggsurv1

## Warning: Removed 5 row(s) containing missing values (geom_path).

## Warning: Removed 4 rows containing missing values (geom_point).

## Warning: Removed 5 row(s) containing missing values (geom_path).

## Warning: Removed 4 rows containing missing values (geom_point).

ggsurv1$plot <- ggsurv1$plot+ 
              ggplot2::annotate("text", 
                                x = 2.5, y = 0.5, # x and y coordinates of the text
                                label = "Adjusted HR (ref: DDLT, HR LDLT 1.20, 95% CI 0.64-2.26)", size = 6)
ggsurv1

## Warning: Removed 5 row(s) containing missing values (geom_path).

## Warning: Removed 4 rows containing missing values (geom_point).

## Warning: Removed 5 row(s) containing missing values (geom_path).

## Warning: Removed 4 rows containing missing values (geom_point).

#Overall
summary(tableby(LDLT ~ Surv(PSURV_years, PCENS), data=UK, times=1:10, surv.stats=c("NeventsSurv", "NriskSurv")))

	0 (N=6498)	1 (N=97)	Total (N=6595)	p value
Surv(PSURV_years, PCENS)				0.686
time = 1	433 (93.2)	8 (91.4)	441 (93.2)
time = 2	579 (90.6)	9 (90.0)	588 (90.6)
time = 3	717 (87.7)	9 (90.0)	726 (87.7)
time = 4	834 (84.7)	10 (88.0)	844 (84.8)
time = 5	917 (82.2)	11 (85.4)	928 (82.3)
time = 6	986 (79.6)	11 (85.4)	997 (79.6)
time = 7	1034 (77.2)	12 (80.0)	1046 (77.3)
time = 8	1074 (74.6)	12 (80.0)	1086 (74.7)
time = 9	1098 (72.4)	13 (66.7)	1111 (72.3)
time = 10	1117 (69.4)	13 (66.7)	1130 (69.4)
time = 1	5531 (93.2)	79 (91.4)	5610 (93.2)
time = 2	4575 (90.6)	65 (90.0)	4640 (90.6)
time = 3	3750 (87.7)	50 (90.0)	3800 (87.7)
time = 4	3059 (84.7)	37 (88.0)	3096 (84.8)
time = 5	2376 (82.2)	27 (85.4)	2403 (82.3)
time = 6	1796 (79.6)	19 (85.4)	1815 (79.6)
time = 7	1319 (77.2)	12 (80.0)	1331 (77.3)
time = 8	945 (74.6)	6 (80.0)	951 (74.7)
time = 9	618 (72.4)	4 (66.7)	622 (72.3)
time = 10	326 (69.4)	2 (66.7)	328 (69.4)

#1 year
summary(tableby(LDLT ~ Surv(PSURV_1year, PCENS_1year), data=UK, times=1:1, surv.stats=c("NeventsSurv", "NriskSurv")))

	0 (N=6498)	1 (N=97)	Total (N=6595)	p value
Surv(PSURV_1year, PCENS_1year)				0.521
time = 1	433 (93.2)	8 (91.4)	441 (93.2)
time = 1	5531 (93.2)	79 (91.4)	5610 (93.2)

#3 year
summary(tableby(LDLT ~ Surv(PSURV_3year, PCENS_3year), data=UK, times=1:3, surv.stats=c("NeventsSurv", "NriskSurv")))

	0 (N=6498)	1 (N=97)	Total (N=6595)	p value
Surv(PSURV_3year, PCENS_3year)				0.692
time = 1	433 (93.2)	8 (91.4)	441 (93.2)
time = 2	579 (90.6)	9 (90.0)	588 (90.6)
time = 3	717 (87.7)	9 (90.0)	726 (87.7)
time = 1	5531 (93.2)	79 (91.4)	5610 (93.2)
time = 2	4575 (90.6)	65 (90.0)	4640 (90.6)
time = 3	3750 (87.7)	50 (90.0)	3800 (87.7)

#5 year
summary(tableby(LDLT ~ Surv(PSURV_5year, PCENS_5year), data=UK, times=1:5, surv.stats=c("NeventsSurv", "NriskSurv")))

	0 (N=6498)	1 (N=97)	Total (N=6595)	p value
Surv(PSURV_5year, PCENS_5year)				0.626
time = 1	433 (93.2)	8 (91.4)	441 (93.2)
time = 2	579 (90.6)	9 (90.0)	588 (90.6)
time = 3	717 (87.7)	9 (90.0)	726 (87.7)
time = 4	834 (84.7)	10 (88.0)	844 (84.8)
time = 5	917 (82.2)	11 (85.4)	928 (82.3)
time = 1	5531 (93.2)	79 (91.4)	5610 (93.2)
time = 2	4575 (90.6)	65 (90.0)	4640 (90.6)
time = 3	3750 (87.7)	50 (90.0)	3800 (87.7)
time = 4	3059 (84.7)	37 (88.0)	3096 (84.8)
time = 5	2376 (82.2)	27 (85.4)	2403 (82.3)

#10 year
summary(tableby(LDLT ~ Surv(PSURV_10year, PCENS_10year), data=UK, times=1:10, surv.stats=c("NeventsSurv", "NriskSurv")))

	0 (N=6498)	1 (N=97)	Total (N=6595)	p value
Surv(PSURV_10year, PCENS_10year)				0.691
time = 1	433 (93.2)	8 (91.4)	441 (93.2)
time = 2	579 (90.6)	9 (90.0)	588 (90.6)
time = 3	717 (87.7)	9 (90.0)	726 (87.7)
time = 4	834 (84.7)	10 (88.0)	844 (84.8)
time = 5	917 (82.2)	11 (85.4)	928 (82.3)
time = 6	986 (79.6)	11 (85.4)	997 (79.6)
time = 7	1034 (77.2)	12 (80.0)	1046 (77.3)
time = 8	1074 (74.6)	12 (80.0)	1086 (74.7)
time = 9	1098 (72.4)	13 (66.7)	1111 (72.3)
time = 10	1117 (69.4)	13 (66.7)	1130 (69.4)
time = 1	5531 (93.2)	79 (91.4)	5610 (93.2)
time = 2	4575 (90.6)	65 (90.0)	4640 (90.6)
time = 3	3750 (87.7)	50 (90.0)	3800 (87.7)
time = 4	3059 (84.7)	37 (88.0)	3096 (84.8)
time = 5	2376 (82.2)	27 (85.4)	2403 (82.3)
time = 6	1796 (79.6)	19 (85.4)	1815 (79.6)
time = 7	1319 (77.2)	12 (80.0)	1331 (77.3)
time = 8	945 (74.6)	6 (80.0)	951 (74.7)
time = 9	618 (72.4)	4 (66.7)	622 (72.3)
time = 10	326 (69.4)	2 (66.7)	328 (69.4)

#Median
survfit(Surv(PSURV_years, PCENS) ~ LDLT, data = UK)

Call: survfit(formula = Surv(PSURV_years, PCENS) ~ LDLT, data = UK)

      n events median 0.95LCL 0.95UCL

LDLT=0 6498 1132 11.9 11.93 NA LDLT=1 97 13 NA 8.13 NA

summary(survfit(Surv(PSURV_years, PCENS)~LDLT, data=UK), times=1:10)

Call: survfit(formula = Surv(PSURV_years, PCENS) ~ LDLT, data = UK)

            LDLT=0

time n.risk n.event survival std.err lower 95% CI upper 95% CI 1 5531 433 0.932 0.00315 0.926 0.938 2 4575 146 0.906 0.00374 0.899 0.913 3 3750 138 0.877 0.00436 0.868 0.885 4 3059 117 0.847 0.00499 0.838 0.857 5 2376 83 0.822 0.00556 0.811 0.833 6 1796 69 0.796 0.00624 0.783 0.808 7 1319 48 0.772 0.00692 0.759 0.786 8 945 40 0.746 0.00783 0.731 0.761 9 618 24 0.724 0.00885 0.706 0.741 10 326 19 0.694 0.01082 0.673 0.716

            LDLT=1

fit4 <- pairwise_survdiff(Surv(PSURV_years, PCENS) ~ LDLT , data = UK)
fit4

Pairwise comparisons using Log-Rank test

data: UK and LDLT

0
1 0.69

P value adjustment method: BH

symnum(fit4$p.value, cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 0.1, 1),
   symbols = c("****", "***", "**", "*", "+", " "),
   abbr.colnames = FALSE, na = "")

0 1
attr(,“legend”) [1] 0 ‘****’ 0.0001 ‘’ 0.001 ’’ 0.01 ’’ 0.05 ‘+’ 0.1 ’ ’ 1

OS overall HCC Figure S2a

HCConly <- bound2 %>% filter(HCC == "1")
fit1 <- survfit(Surv(PSURV_years, PCENS) ~ COUNTRY, data = HCConly)

ggsurv1 <- ggsurvplot(
           fit1,               
           data = HCConly,          
           risk.table = TRUE,     
           pval = TRUE,
           pval.coord = c(0, 0.8),
           pval.size = 6,
           conf.int = F,         
           xlim = c(0,5),
           ylim = c(0.7, 1.00),
           xlab = "Years from LT",
           palette = "jama", #dont use red and green
           break.time.by = 1,  
           title = "",
      # subtitle = "with 95% confidence intervals",
               ggtheme = theme_test(base_size=28, base_family = "Helvetica"),
         # risk.table.y.text.col = F,
        fontsize=8,
         risk.table.height = 0.25,
        #  risk.table.y.text = T,
        #  surv.median.line = "hv", 
           legend.title= "",
      tables.x.text = FALSE,
          legend.labs =
           c("US", "Canada", "UK")
        )

## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.

ggsurv1$plot <- ggsurv1$plot+ 
              ggplot2::annotate("text", 
                                x = 2.5, y = 0.7, # x and y coordinates of the text
                                label = "Adjusted HR (ref: US, Canada HR 0.20, 95% CI 0.02-1.66, UK HR 0.38, 95% CI 0.05-2.96)", size = 6)
ggsurv1

## Warning: Removed 84 row(s) containing missing values (geom_path).

## Warning: Removed 71 rows containing missing values (geom_point).

## Warning: Removed 84 row(s) containing missing values (geom_path).

## Warning: Removed 71 rows containing missing values (geom_point).

OS overall non-HCC Figure S2b

nonHCC <- bound2 %>% filter(HCC == "0")
fit1 <- survfit(Surv(PSURV_years, PCENS) ~ COUNTRY, data = nonHCC)

ggsurv1 <- ggsurvplot(
           fit1,               
           data = nonHCC,          
           risk.table = TRUE,     
           pval = TRUE,
           pval.coord = c(0, 0.8),
           pval.size = 6,
           conf.int = F,         
           xlim = c(0,5),
           ylim = c(0.7, 1.00),
           xlab = "Years from LT",
           palette = "jama",
           break.time.by = 1,  
           title = "",
      # subtitle = "with 95% confidence intervals",
               ggtheme = theme_test(base_size=28, base_family = "Helvetica"),
         # risk.table.y.text.col = F,
         fontsize=8,
         risk.table.height = 0.25,
        #  risk.table.y.text = T,
        #  surv.median.line = "hv", 
           legend.title= "",
      tables.x.text = FALSE,
          legend.labs =
           c("US", "Canada", "UK")
        )

## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.

ggsurv1

## Warning: Removed 15 row(s) containing missing values (geom_path).

## Warning: Removed 14 rows containing missing values (geom_point).

## Warning: Removed 15 row(s) containing missing values (geom_path).

## Warning: Removed 14 rows containing missing values (geom_point).

ggsurv1$plot <- ggsurv1$plot+ 
              ggplot2::annotate("text", 
                                x = 2.5, y = 0.7, # x and y coordinates of the text
                                label = "Adjusted HR (ref: US, Canada HR 0.51, 95% CI 0.22-1.18, UK HR 1.36, 95% CI 0.69-2.67)", size = 6)
ggsurv1

## Warning: Removed 15 row(s) containing missing values (geom_path).

## Warning: Removed 14 rows containing missing values (geom_point).

## Warning: Removed 15 row(s) containing missing values (geom_path).

## Warning: Removed 14 rows containing missing values (geom_point).

confint.coxme <- function(object, parm=NULL, level=0.95, ..., more=FALSE){
  if(!is.null(parm)) warning("[confint.coxme] argument 'parm' doesn't do anything for this method")
  if(level != 0.95) warning("[confint.coxme] 'level' will be 0.95 regardless of what argument you give it. Ha!")
  co <- object$coef
  se <- sqrt(diag(stats::vcov(object)))
  m <- matrix(c(co - 2*se, co + 2*se), ncol=2)
  colnames(m) <- c("2.5 %", "97.5 %")
  rownames(m) <- names(co)
  if(more){
    p <- 2*stats::pnorm(abs(co/se), lower.tail=F)
    m <- cbind(m, co, p)
    rownames(m)[3:4] <- c("coef", "p")
  }
  return (m)
}

US vs. UK vs. CAN Forest model patient survival

bound4 <- bound2 %>% filter(COUNTRY == "US" | COUNTRY == "UK" | COUNTRY == "CAN")

#Uni
cox1 <- coxph(Surv(PSURV,PCENS)~ COUNTRY, data=bound4)
forest_model(cox1)

#Just recipient
cox1 <- coxme(formula= Surv(PSURV_years,PCENS)~COUNTRY+
              #RSEX+
              BMI+  
              RBILIRUBIN+
              RINR+
              RCREAT+
              RREN_SUP+
              RAGE+
              WAITLIST_TIME+
              TX_YR+
              HCV+
              HCC+
              NASH+
              ALD+
              (1|TRANSPLANT_UNIT),
        data=bound4)
cox1

Cox mixed-effects model fit by maximum likelihood Data: bound4 events, n = 414, 2692 (394 observations deleted due to missingness) Iterations= 11 80 NULL Integrated Fitted Log-likelihood -3040.108 -3003.128 -2996.95

              Chisq    df                p   AIC    BIC

Integrated loglik 73.96 15.00 0.00000000087150 43.96 -16.43 Penalized loglik 86.32 19.44 0.00000000020941 47.44 -30.81

Model: Surv(PSURV_years, PCENS) ~ COUNTRY + BMI + RBILIRUBIN + RINR + RCREAT + RREN_SUP + RAGE + WAITLIST_TIME + TX_YR + HCV + HCC + NASH + ALD + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z p COUNTRYCAN -0.5311240589 0.5879437 0.2375154734 -2.24 0.025000 COUNTRYUK 0.2855994130 1.3305593 0.3080447373 0.93 0.350000 BMI -0.0087710805 0.9912673 0.0103955919 -0.84 0.400000 RBILIRUBIN -0.0090471908 0.9909936 0.0122679672 -0.74 0.460000 RINR -0.0476529127 0.9534647 0.0963508969 -0.49 0.620000 RCREAT 0.2525526779 1.2873073 0.0741830872 3.40 0.000660 RREN_SUPPre-tx support 0.8457528619 2.3297311 0.4370902360 1.93 0.053000 RAGE 0.0189438183 1.0191244 0.0048279739 3.92 0.000087 WAITLIST_TIME -0.0000445156 0.9999555 0.0000972877 -0.46 0.650000 TX_YR -0.0348418882 0.9657581 0.0191004486 -1.82 0.068000 HCV1 -0.6577760333 0.5180021 0.3290465504 -2.00 0.046000 HCC1 0.3612829760 1.4351695 0.1244709393 2.90 0.003700 NASH1 -0.0460932088 0.9549529 0.1574176358 -0.29 0.770000 ALD1 -0.1821803103 0.8334510 0.1841702976 -0.99 0.320000

Random effects Group Variable Std Dev Variance
TRANSPLANT_UNIT Intercept 0.13802843 0.01905185

confint.coxme(cox1)

                          2.5 %        97.5 %

COUNTRYCAN -1.006155006 -0.0560931122 COUNTRYUK -0.330490062 0.9016888875 BMI -0.029562264 0.0120201032 RBILIRUBIN -0.033583125 0.0154887436 RINR -0.240354707 0.1450488811 RCREAT 0.104186503 0.4009188524 RREN_SUPPre-tx support -0.028427610 1.7199333340 RAGE 0.009287871 0.0285997661 WAITLIST_TIME -0.000239091 0.0001500598 TX_YR -0.073042785 0.0033590090 HCV1 -1.315869134 0.0003170674 HCC1 0.112341097 0.6102248546 NASH1 -0.360928480 0.2687420627 ALD1 -0.550520905 0.1861602848

#Just donor
cox1 <- coxme(formula= Surv(PSURV_years,PCENS)~COUNTRY+
              CIT+
              DAGE+
              DSEX+
              TX_YR+
              (1|TRANSPLANT_UNIT),
        data=bound4)
cox1

Cox mixed-effects model fit by maximum likelihood Data: bound4 events, n = 389, 2756 (330 observations deleted due to missingness) Iterations= 11 58 NULL Integrated Fitted Log-likelihood -2856.467 -2836.57 -2832.481

              Chisq  df             p   AIC   BIC

Integrated loglik 39.79 7.0 0.00000137840 25.79 -1.95 Penalized loglik 47.97 9.7 0.00000048396 28.57 -9.87

Model: Surv(PSURV_years, PCENS) ~ COUNTRY + CIT + DAGE + DSEX + TX_YR + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z p COUNTRYCAN -0.81837773001 0.4411467 0.2241123247 -3.65 0.000260 COUNTRYUK -0.00073967936 0.9992606 0.3006234626 0.00 1.000000 CIT 0.00008032782 1.0000803 0.0002061973 0.39 0.700000 DAGE 0.01993839259 1.0201385 0.0048941671 4.07 0.000046 DSEXMale 0.08047482060 1.0838016 0.1031549993 0.78 0.440000 TX_YR -0.02327597647 0.9769928 0.0191499160 -1.22 0.220000

Random effects Group Variable Std Dev Variance
TRANSPLANT_UNIT Intercept 0.11222883 0.01259531

confint.coxme(cox1)

               2.5 %        97.5 %

COUNTRYCAN -1.2666023795 -0.3701530805 COUNTRYUK -0.6019866045 0.6005072458 CIT -0.0003320668 0.0004927224 DAGE 0.0101500583 0.0297267268 DSEXMale -0.1258351780 0.2867848192 TX_YR -0.0615758085 0.0150238556

cox2 <- coxme(formula= Surv(PSURV_years,PCENS)~COUNTRY+
              #RSEX+
              BMI+  
              RBILIRUBIN+
              RINR+
              RCREAT+
              RREN_SUP+
              RAGE+
              WAITLIST_TIME+
              TX_YR+
              HCV+
              HCC+
              NASH+
              ALD+
              CIT+
              DAGE+
              DSEX+
               (1|TRANSPLANT_UNIT),
        data=bound4)
cox2

Cox mixed-effects model fit by maximum likelihood Data: bound4 events, n = 370, 2454 (632 observations deleted due to missingness) Iterations= 15 108 NULL Integrated Fitted Log-likelihood -2672.585 -2630.344 -2625.967

              Chisq    df                 p   AIC    BIC

Integrated loglik 84.48 18.00 0.000000000139260 48.48 -21.96 Penalized loglik 93.24 20.94 0.000000000042409 51.35 -30.60

Model: Surv(PSURV_years, PCENS) ~ COUNTRY + BMI + RBILIRUBIN + RINR + RCREAT + RREN_SUP + RAGE + WAITLIST_TIME + TX_YR + HCV + HCC + NASH + ALD + CIT + DAGE + DSEX + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z p COUNTRYCAN -0.63572243800 0.5295528 0.38285523502 -1.66 0.097000 COUNTRYUK 0.29150188194 1.3384362 0.31812715216 0.92 0.360000 BMI -0.03378677979 0.9667776 0.01174515328 -2.88 0.004000 RBILIRUBIN -0.00500563594 0.9950069 0.01298705057 -0.39 0.700000 RINR -0.04096760903 0.9598602 0.09891411935 -0.41 0.680000 RCREAT 0.28150105342 1.3251174 0.08429903733 3.34 0.000840 RREN_SUPPre-tx support 0.60753251684 1.8358958 0.47870798106 1.27 0.200000 RAGE 0.02163205949 1.0218677 0.00504329523 4.29 0.000018 WAITLIST_TIME -0.00003437378 0.9999656 0.00009670316 -0.36 0.720000 TX_YR -0.03111807712 0.9693611 0.02021414227 -1.54 0.120000 HCV1 -1.18034770023 0.3071719 0.51036105217 -2.31 0.021000 HCC1 0.30786328749 1.3605150 0.13355522904 2.31 0.021000 NASH1 0.04476354781 1.0457806 0.16283139068 0.27 0.780000 ALD1 -0.07171426584 0.9307968 0.18612856088 -0.39 0.700000 CIT 0.00006667081 1.0000667 0.00020777368 0.32 0.750000 DAGE 0.01666608413 1.0168057 0.00514656037 3.24 0.001200 DSEXMale 0.09010499534 1.0942892 0.10635731764 0.85 0.400000

Random effects Group Variable Std Dev Variance
TRANSPLANT_UNIT Intercept 0.12042989 0.01450336

confint.coxme(cox2)

                           2.5 %        97.5 %

COUNTRYCAN -1.4014329080 0.1299880320 COUNTRYUK -0.3447524224 0.9277561863 BMI -0.0572770864 -0.0102964732 RBILIRUBIN -0.0309797371 0.0209684652 RINR -0.2387958477 0.1568606297 RCREAT 0.1129029788 0.4500991281 RREN_SUPPre-tx support -0.3498834453 1.5649484790 RAGE 0.0115454690 0.0317186499 WAITLIST_TIME -0.0002277801 0.0001590325 TX_YR -0.0715463617 0.0093102074 HCV1 -2.2010698046 -0.1596255959 HCC1 0.0407528294 0.5749737456 NASH1 -0.2808992336 0.3704263292 ALD1 -0.4439713876 0.3005428559 CIT -0.0003488765 0.0004822182 DAGE 0.0063729634 0.0269592049 DSEXMale -0.1226096399 0.3028196306

US vs. UK vs. CAN Forest model graft survival

bound4 <- bound2 %>% filter(COUNTRY == "US" | COUNTRY == "UK" | COUNTRY == "CAN")

#Uni
cox1 <- coxph(Surv(GSURV_years, GCENS)~ COUNTRY, data=bound4)
forest_model(cox1)

#Just recipient
cox1 <- coxme(formula= Surv(GSURV_years, GCENS)~COUNTRY+
              #RSEX+
              BMI+  
              RBILIRUBIN+
              RINR+
              RCREAT+
              RREN_SUP+
              RAGE+
              WAITLIST_TIME+
              TX_YR+
              HCV+
              HCC+
              NASH+
              ALD+
              (1|TRANSPLANT_UNIT),
        data=bound4)
cox1

Cox mixed-effects model fit by maximum likelihood Data: bound4 events, n = 578, 2692 (394 observations deleted due to missingness) Iterations= 5 33 NULL Integrated Fitted Log-likelihood -4293.493 -4266.617 -4260.902

              Chisq    df             p   AIC    BIC

Integrated loglik 53.75 15.00 0.00000288840 23.75 -41.64 Penalized loglik 65.18 19.11 0.00000061248 26.97 -56.34

Model: Surv(GSURV_years, GCENS) ~ COUNTRY + BMI + RBILIRUBIN + RINR + RCREAT + RREN_SUP + RAGE + WAITLIST_TIME + TX_YR + HCV + HCC + NASH + ALD + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z p COUNTRYCAN -0.610558564149 0.5430475 0.2094216039 -2.92 0.00360 COUNTRYUK -0.126940611306 0.8807860 0.3015188443 -0.42 0.67000 BMI -0.033234897047 0.9673113 0.0094756215 -3.51 0.00045 RBILIRUBIN -0.000651288247 0.9993489 0.0088870342 -0.07 0.94000 RINR 0.007262250223 1.0072887 0.0630546708 0.12 0.91000 RCREAT 0.219535233446 1.2454977 0.0716568453 3.06 0.00220 RREN_SUPPre-tx support 0.273768216684 1.3149100 0.4460513851 0.61 0.54000 RAGE -0.000176919462 0.9998231 0.0036667537 -0.05 0.96000 WAITLIST_TIME 0.000001032512 1.0000010 0.0000761251 0.01 0.99000 TX_YR -0.030819085407 0.9696510 0.0154761589 -1.99 0.04600 HCV1 -0.819829727009 0.4405067 0.3249926519 -2.52 0.01200 HCC1 0.280171824527 1.3233572 0.1116098632 2.51 0.01200 NASH1 0.098851075646 1.1039019 0.1365813115 0.72 0.47000 ALD1 -0.073951936862 0.9287163 0.1514441783 -0.49 0.63000

Random effects Group Variable Std Dev Variance
TRANSPLANT_UNIT Intercept 0.11106823 0.01233615

confint.coxme(cox1)

                           2.5 %        97.5 %

COUNTRYCAN -1.0294017720 -0.1917153563 COUNTRYUK -0.7299782999 0.4760970773 BMI -0.0521861401 -0.0142836540 RBILIRUBIN -0.0184253567 0.0171227802 RINR -0.1188470915 0.1333715919 RCREAT 0.0762215429 0.3628489240 RREN_SUPPre-tx support -0.6183345536 1.1658709869 RAGE -0.0075104268 0.0071565878 WAITLIST_TIME -0.0001512177 0.0001532827 TX_YR -0.0617714031 0.0001332323 HCV1 -1.4698150307 -0.1698444233 HCC1 0.0569520982 0.5033915509 NASH1 -0.1743115474 0.3720136987 ALD1 -0.3768402934 0.2289364197

#Just donor
cox1 <- coxme(formula= Surv(GSURV_years, GCENS)~COUNTRY+
              CIT+
              DAGE+
              DSEX+
              TX_YR+
              (1|TRANSPLANT_UNIT),
        data=bound4)
cox1

Cox mixed-effects model fit by maximum likelihood Data: bound4 events, n = 550, 2756 (330 observations deleted due to missingness) Iterations= 5 27 NULL Integrated Fitted Log-likelihood -4090.247 -4063.92 -4060.432

              Chisq   df               p   AIC  BIC

Integrated loglik 52.65 7.00 0.0000000043389 38.65 8.48 Penalized loglik 59.63 9.18 0.0000000018934 41.27 1.72

Model: Surv(GSURV_years, GCENS) ~ COUNTRY + CIT + DAGE + DSEX + TX_YR + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z p COUNTRYCAN -0.72102756584 0.4862523 0.1830105943 -3.94 0.0000820 COUNTRYUK -0.29139845920 0.7472179 0.2966431756 -0.98 0.3300000 CIT -0.00002391839 0.9999761 0.0001965816 -0.12 0.9000000 DAGE 0.02006867245 1.0202714 0.0041112230 4.88 0.0000011 DSEXMale 0.00791046142 1.0079418 0.0867379149 0.09 0.9300000 TX_YR -0.03139422217 0.9690935 0.0154944159 -2.03 0.0430000

Random effects Group Variable Std Dev Variance
TRANSPLANT_UNIT Intercept 0.086267299 0.007442047

confint.coxme(cox1)

               2.5 %        97.5 %

COUNTRYCAN -1.0870487544 -0.3550063773 COUNTRYUK -0.8846848104 0.3018878919 CIT -0.0004170816 0.0003692448 DAGE 0.0118462265 0.0282911184 DSEXMale -0.1655653683 0.1813862912 TX_YR -0.0623830540 -0.0004053904

cox2 <- coxme(formula= Surv(GSURV_years, GCENS)~COUNTRY+
              #RSEX+
              BMI+  
              RBILIRUBIN+
              RINR+
              RCREAT+
              RREN_SUP+
              RAGE+
              WAITLIST_TIME+
              TX_YR+
              HCV+
              HCC+
              NASH+
              ALD+
              CIT+
              DAGE+
              DSEX+
               (1|TRANSPLANT_UNIT),
        data=bound4)
cox2

Cox mixed-effects model fit by maximum likelihood Data: bound4 events, n = 522, 2454 (632 observations deleted due to missingness) Iterations= 5 32 NULL Integrated Fitted Log-likelihood -3824.806 -3791.884 -3788.624

              Chisq    df              p   AIC    BIC

Integrated loglik 65.84 18.00 0.000000225150 29.84 -46.79 Penalized loglik 72.36 19.99 0.000000073919 32.39 -52.71

Model: Surv(GSURV_years, GCENS) ~ COUNTRY + BMI + RBILIRUBIN + RINR + RCREAT + RREN_SUP + RAGE + WAITLIST_TIME + TX_YR + HCV + HCC + NASH + ALD + CIT + DAGE + DSEX + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z p COUNTRYCAN -0.607031225699 0.5449664 0.30955812759 -1.96 0.05000 COUNTRYUK -0.096185415670 0.9082956 0.31216577522 -0.31 0.76000 BMI -0.035293132541 0.9653224 0.00995385102 -3.55 0.00039 RBILIRUBIN 0.005831617010 1.0058487 0.00923977530 0.63 0.53000 RINR -0.000788312547 0.9992120 0.06828058260 -0.01 0.99000 RCREAT 0.235429885844 1.2654527 0.08122336231 2.90 0.00370 RREN_SUPPre-tx support 0.008258619433 1.0082928 0.50200452296 0.02 0.99000 RAGE 0.002632655953 1.0026361 0.00388534503 0.68 0.50000 WAITLIST_TIME -0.000007797856 0.9999922 0.00007801948 -0.10 0.92000 TX_YR -0.036091403477 0.9645521 0.01619140184 -2.23 0.02600 HCV1 -1.161818804516 0.3129165 0.45436992157 -2.56 0.01100 HCC1 0.258328576372 1.2947642 0.11740904739 2.20 0.02800 NASH1 0.104915180452 1.1106164 0.14075116123 0.75 0.46000 ALD1 -0.038070621089 0.9626450 0.15562983286 -0.24 0.81000 CIT -0.000007497090 0.9999925 0.00019703100 -0.04 0.97000 DAGE 0.018274758472 1.0184428 0.00427990577 4.27 0.00002 DSEXMale 0.019382431009 1.0195715 0.08933865860 0.22 0.83000

Random effects Group Variable Std Dev Variance
TRANSPLANT_UNIT Intercept 0.085642010 0.007334554

confint.coxme(cox2)

                           2.5 %        97.5 %

COUNTRYCAN -1.2261474809 0.0120850295 COUNTRYUK -0.7205169661 0.5281461348 BMI -0.0552008346 -0.0153854305 RBILIRUBIN -0.0126479336 0.0243111676 RINR -0.1373494777 0.1357728526 RCREAT 0.0729831612 0.3978766105 RREN_SUPPre-tx support -0.9957504265 1.0122676654 RAGE -0.0051380341 0.0104033460 WAITLIST_TIME -0.0001638368 0.0001482411 TX_YR -0.0684742072 -0.0037085998 HCV1 -2.0705586476 -0.2530789614 HCC1 0.0235104816 0.4931466712 NASH1 -0.1765871420 0.3864175029 ALD1 -0.3493302868 0.2731890446 CIT -0.0004015591 0.0003865649 DAGE 0.0097149469 0.0268345700 DSEXMale -0.1592948862 0.1980597482

#Uni
cox1 <- coxph(Surv(PSURV,PCENS)~ LDLT, data=US)
forest_model(cox1)

## Resized limits to included dashed line in forest panel

#Just recipient
cox1 <- coxme(formula= Surv(PSURV_years,PCENS)~LDLT+
             # RSEX+
              BMI+  
              RBILIRUBIN+
              RINR+
              RCREAT+
              RREN_SUP+
              RAGE+
              WAITLIST_TIME+
              TX_YR+
              HCV+
              HCC+
              NASH+
              ALD+
              (1|TRANSPLANT_UNIT),
        data=US)
cox1

Cox mixed-effects model fit by maximum likelihood Data: US events, n = 13054, 60055 (130 observations deleted due to missingness) Iterations= 18 94 NULL Integrated Fitted Log-likelihood -136255.1 -135580.3 -135461.9

                Chisq    df p     AIC     BIC

Integrated loglik 1349.68 14.00 0 1321.68 1217.00 Penalized loglik 1586.40 93.65 0 1399.10 698.91

Model: Surv(PSURV_years, PCENS) ~ LDLT + BMI + RBILIRUBIN + RINR + RCREAT + RREN_SUP + RAGE + WAITLIST_TIME + TX_YR + HCV + HCC + NASH + ALD + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z p LDLT1 -0.138396072553 0.8707537 0.05482507482 -2.52 0.01200 BMI -0.005333791129 0.9946804 0.00162002063 -3.29 0.00099 RBILIRUBIN 0.002358194938 1.0023610 0.00093278365 2.53 0.01100 RINR -0.002658872265 0.9973447 0.00741570882 -0.36 0.72000 RCREAT 0.079546683215 1.0827961 0.00808780178 9.84 0.00000 RREN_SUPPre-tx support 0.344891007973 1.4118360 0.03230276716 10.68 0.00000 RAGE 0.020150945955 1.0203553 0.00103200067 19.53 0.00000 WAITLIST_TIME 0.000006135301 1.0000061 0.00001796821 0.34 0.73000 TX_YR -0.049904919977 0.9513199 0.00337978386 -14.77 0.00000 HCV1 -0.001301056103 0.9986998 0.03952833720 -0.03 0.97000 HCC1 0.201305676124 1.2229986 0.02155153884 9.34 0.00000 NASH1 -0.027697407178 0.9726826 0.02787684515 -0.99 0.32000 ALD1 -0.080163957682 0.9229650 0.02908264164 -2.76 0.00580

Random effects Group Variable Std Dev Variance
TRANSPLANT_UNIT Intercept 0.17625600 0.03106618

confint.coxme(cox1)

                            2.5 %         97.5 %

LDLT1 -0.24804622219 -0.02874592292 BMI -0.00857383239 -0.00209374987 RBILIRUBIN 0.00049262763 0.00422376224 RINR -0.01749028990 0.01217254537 RCREAT 0.06337107966 0.09572228677 RREN_SUPPre-tx support 0.28028547365 0.40949654230 RAGE 0.01808694461 0.02221494730 WAITLIST_TIME -0.00002980112 0.00004207172 TX_YR -0.05666448769 -0.04314535226 HCV1 -0.08035773050 0.07775561829 HCC1 0.15820259845 0.24440875380 NASH1 -0.08345109748 0.02805628312 ALD1 -0.13832924095 -0.02199867441

#Just donor
cox1 <- coxme(formula= Surv(PSURV_years,PCENS)~LDLT+
              CIT+
              DAGE+
              DSEX+
              TX_YR+
              (1|TRANSPLANT_UNIT),
        data=US)
cox1

Cox mixed-effects model fit by maximum likelihood Data: US events, n = 12924, 59695 (490 observations deleted due to missingness) Iterations= 17 72 NULL Integrated Fitted Log-likelihood -134809.3 -134480.3 -134365.8

               Chisq    df p    AIC    BIC

Integrated loglik 658.02 6.00 0 646.02 601.22 Penalized loglik 887.02 84.18 0 718.67 90.14

Model: Surv(PSURV_years, PCENS) ~ LDLT + CIT + DAGE + DSEX + TX_YR + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z p LDLT1 -0.1428864758 0.8668525 0.05746830508 -2.49 0.0130000000 CIT 0.0002993341 1.0002994 0.00005066372 5.91 0.0000000035 DAGE 0.0079183288 1.0079498 0.00054805557 14.45 0.0000000000 DSEXMale 0.0333538423 1.0339163 0.01817504103 1.84 0.0660000000 TX_YR -0.0408393881 0.9599833 0.00334509653 -12.21 0.0000000000

Random effects Group Variable Std Dev Variance
TRANSPLANT_UNIT Intercept 0.17015352 0.02895222

confint.coxme(cox1)

             2.5 %        97.5 %

LDLT1 -0.2578230860 -0.0279498657 CIT 0.0001980067 0.0004006615 DAGE 0.0068222176 0.0090144399 DSEXMale -0.0029962398 0.0697039243 TX_YR -0.0475295812 -0.0341491951

#donor and recipient
cox2 <- coxme(formula= Surv(PSURV_years,PCENS)~LDLT+
              #RSEX+
              BMI+  
              RBILIRUBIN+
              RINR+
              RCREAT+
              RREN_SUP+
              RAGE+
              WAITLIST_TIME+
              TX_YR+
              HCV+
              HCC+
              NASH+
              ALD+
              CIT+
              DAGE+
              DSEX+
               (1|TRANSPLANT_UNIT),
        data=US)
cox2

Cox mixed-effects model fit by maximum likelihood Data: US events, n = 12893, 59567 (618 observations deleted due to missingness) Iterations= 19 99 NULL Integrated Fitted Log-likelihood -134463.3 -133688.9 -133571.1

                Chisq    df p     AIC     BIC

Integrated loglik 1548.76 17.00 0 1514.76 1387.86 Penalized loglik 1784.42 96.36 0 1591.70 872.40

Model: Surv(PSURV_years, PCENS) ~ LDLT + BMI + RBILIRUBIN + RINR + RCREAT + RREN_SUP + RAGE + WAITLIST_TIME + TX_YR + HCV + HCC + NASH + ALD + CIT + DAGE + DSEX + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z LDLT1 0.0310905119085 1.0315789 0.05834312860 0.53 BMI -0.0060951129791 0.9939234 0.00163570893 -3.73 RBILIRUBIN 0.0029983182905 1.0030028 0.00093552912 3.20 RINR -0.0000705908230 0.9999294 0.00723002548 -0.01 RCREAT 0.0805673550288 1.0839019 0.00809120861 9.96 RREN_SUPPre-tx support 0.3654363945530 1.4411428 0.03249890432 11.24 RAGE 0.0192502364733 1.0194367 0.00104292529 18.46 WAITLIST_TIME 0.0000007959007 1.0000008 0.00001809947 0.04 TX_YR -0.0471220675239 0.9539709 0.00342447487 -13.76 HCV1 -0.0035821793134 0.9964242 0.03976106539 -0.09 HCC1 0.2069491588403 1.2299200 0.02167981000 9.55 NASH1 -0.0405861375053 0.9602264 0.02804873689 -1.45 ALD1 -0.0960670891460 0.9084031 0.02925449134 -3.28 CIT 0.0003296642817 1.0003297 0.00005011103 6.58 DAGE 0.0076736704493 1.0077032 0.00055385747 13.85 DSEXMale 0.0372574774604 1.0379602 0.01823068689 2.04 p LDLT1 0.590000000000 BMI 0.000190000000 RBILIRUBIN 0.001400000000 RINR 0.990000000000 RCREAT 0.000000000000 RREN_SUPPre-tx support 0.000000000000 RAGE 0.000000000000 WAITLIST_TIME 0.960000000000 TX_YR 0.000000000000 HCV1 0.930000000000 HCC1 0.000000000000 NASH1 0.150000000000 ALD1 0.001000000000 CIT 0.000000000047 DAGE 0.000000000000 DSEXMale 0.041000000000

Random effects Group Variable Std Dev Variance TRANSPLANT_UNIT Intercept 0.1764769 0.0311441

confint.coxme(cox2)

                            2.5 %         97.5 %

LDLT1 -0.08559574530 0.14777676912 BMI -0.00936653084 -0.00282369512 RBILIRUBIN 0.00112726004 0.00486937654 RINR -0.01453064178 0.01438946014 RCREAT 0.06438493781 0.09674977224 RREN_SUPPre-tx support 0.30043858591 0.43043420319 RAGE 0.01716438588 0.02133608706 WAITLIST_TIME -0.00003540303 0.00003699483 TX_YR -0.05397101726 -0.04027311779 HCV1 -0.08310431010 0.07593995147 HCC1 0.16358953883 0.25030877885 NASH1 -0.09668361129 0.01551133628 ALD1 -0.15457607182 -0.03755810648 CIT 0.00022944223 0.00042988633 DAGE 0.00656595551 0.00878138539 DSEXMale 0.00079610367 0.07371885125

#Uni
cox1 <- coxph(Surv(PSURV,PCENS)~ LDLT, data=CAN)
forest_model(cox1)

## Resized limits to included dashed line in forest panel

#Just recipient
cox1 <- coxme(formula= Surv(PSURV_years,PCENS)~LDLT+
             # RSEX+
              BMI+  
              RBILIRUBIN+
              RINR+
              RCREAT+
              RREN_SUP+
              RAGE+
              WAITLIST_TIME+
              TX_YR+
              HCV+
              HCC+
              NASH+
              ALD+
              (1|TRANSPLANT_UNIT),
        data=CAN)
cox1

Cox mixed-effects model fit by maximum likelihood Data: CAN events, n = 275, 1562 (2064 observations deleted due to missingness) Iterations= 8 42 NULL Integrated Fitted Log-likelihood -1911.173 -1894.755 -1894.754

              Chisq df         p  AIC    BIC

Integrated loglik 32.84 14 0.0030410 4.84 -45.80 Penalized loglik 32.84 13 0.0018026 6.84 -40.19

Model: Surv(PSURV_years, PCENS) ~ LDLT + BMI + RBILIRUBIN + RINR + RCREAT + RREN_SUP + RAGE + WAITLIST_TIME + TX_YR + HCV + HCC + NASH + ALD + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z p LDLT1 -0.4059059463 0.6663728 0.2139677437 -1.90 0.0580 BMI 0.0070975297 1.0071228 0.0047597144 1.49 0.1400 RBILIRUBIN -0.0009452539 0.9990552 0.0058360665 -0.16 0.8700 RINR 0.0830430941 1.0865886 0.0662524387 1.25 0.2100 RCREAT 0.0900951266 1.0942784 0.0623984713 1.44 0.1500 RREN_SUPPre-tx support 0.0521201370 1.0535023 0.6229937905 0.08 0.9300 RAGE 0.0213351425 1.0215644 0.0070187279 3.04 0.0024 WAITLIST_TIME -0.0000195753 0.9999804 0.0001664136 -0.12 0.9100 TX_YR -0.0053806576 0.9946338 0.0240358638 -0.22 0.8200 HCV1 0.2394025201 1.2704898 0.1347490888 1.78 0.0760 HCC1 0.1125268684 1.1191023 0.1453354977 0.77 0.4400 NASH1 0.1353453565 1.1449321 0.2303036131 0.59 0.5600 ALD1 -0.0874127530 0.9162988 0.1560258829 -0.56 0.5800

Random effects Group Variable Std Dev Variance
TRANSPLANT_UNIT Intercept 0.00397686599 0.00001581546

confint.coxme(cox1)

                           2.5 %       97.5 %

LDLT1 -0.8338414337 0.0220295411 BMI -0.0024218990 0.0166169585 RBILIRUBIN -0.0126173868 0.0107268790 RINR -0.0494617832 0.2155479715 RCREAT -0.0347018161 0.2148920693 RREN_SUPPre-tx support -1.1938674439 1.2981077180 RAGE 0.0072976867 0.0353725982 WAITLIST_TIME -0.0003524025 0.0003132519 TX_YR -0.0534523852 0.0426910700 HCV1 -0.0300956576 0.5089006977 HCC1 -0.1781441270 0.4031978638 NASH1 -0.3252618697 0.5959525826 ALD1 -0.3994645188 0.2246390129

#Just donor
cox1 <- coxme(formula= Surv(PSURV_years,PCENS)~LDLT+
              CIT+
              DAGE+
              DSEX+
              TX_YR+
              (1|TRANSPLANT_UNIT),
        data=CAN)
cox1

Cox mixed-effects model fit by maximum likelihood Data: CAN events, n = 387, 2717 (909 observations deleted due to missingness) Iterations= 7 37 NULL Integrated Fitted Log-likelihood -2896.113 -2877.986 -2877.181

              Chisq   df             p   AIC  BIC

Integrated loglik 36.25 6.00 0.00000246110 24.25 0.50 Penalized loglik 37.86 5.61 0.00000079486 26.63 4.41

Model: Surv(PSURV_years, PCENS) ~ LDLT + CIT + DAGE + DSEX + TX_YR + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z p LDLT1 -0.5285835670 0.5894393 0.2274256290 -2.32 0.0200 CIT 0.0005932012 1.0005934 0.0002719127 2.18 0.0290 DAGE 0.0088863681 1.0089260 0.0030151317 2.95 0.0032 DSEXMale 0.0702133296 1.0727370 0.1138151879 0.62 0.5400 TX_YR -0.0192695022 0.9809150 0.0200001840 -0.96 0.3400

Random effects Group Variable Std Dev Variance
TRANSPLANT_UNIT Intercept 0.095737942 0.009165753

confint.coxme(cox1)

              2.5 %       97.5 %

LDLT1 -0.98343482505 -0.073732309 CIT 0.00004937588 0.001137027 DAGE 0.00285610479 0.014916631 DSEXMale -0.15741704618 0.297843705 TX_YR -0.05926987022 0.020730866

#donor and recipient
cox2 <- coxme(formula= Surv(PSURV_years,PCENS)~LDLT+
             # RSEX+
              BMI+  
              RBILIRUBIN+
              RINR+
              RCREAT+
              RREN_SUP+
              RAGE+
              WAITLIST_TIME+
              TX_YR+
              HCV+
              HCC+
              NASH+
              ALD+
              CIT+
              DAGE+
              DSEX+
               (1|TRANSPLANT_UNIT),
        data=CAN)
cox2

Cox mixed-effects model fit by maximum likelihood Data: CAN events, n = 209, 1174 (2452 observations deleted due to missingness) Iterations= 5 27 NULL Integrated Fitted Log-likelihood -1382.751 -1363.347 -1363.345

              Chisq df         p  AIC    BIC

Integrated loglik 38.81 17 0.0018995 4.81 -52.01 Penalized loglik 38.81 16 0.0011587 6.81 -46.68

Model: Surv(PSURV_years, PCENS) ~ LDLT + BMI + RBILIRUBIN + RINR + RCREAT + RREN_SUP + RAGE + WAITLIST_TIME + TX_YR + HCV + HCC + NASH + ALD + CIT + DAGE + DSEX + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z p LDLT1 -0.42375444431 0.6545846 0.3838033773 -1.10 0.2700 BMI 0.00256551883 1.0025688 0.0075680658 0.34 0.7300 RBILIRUBIN 0.00517942937 1.0051929 0.0062152900 0.83 0.4000 RINR 0.06956679059 1.0720437 0.0746042913 0.93 0.3500 RCREAT 0.03016805276 1.0306277 0.0736548232 0.41 0.6800 RREN_SUPPre-tx support -0.69335897371 0.4998941 1.0321275595 -0.67 0.5000 RAGE 0.02465131144 1.0249577 0.0082866567 2.97 0.0029 WAITLIST_TIME 0.00006299711 1.0000630 0.0001813164 0.35 0.7300 TX_YR 0.01008658292 1.0101376 0.0277848762 0.36 0.7200 HCV1 0.25178213497 1.2863158 0.1547902429 1.63 0.1000 HCC1 0.05001748240 1.0512895 0.1659043526 0.30 0.7600 NASH1 0.20608716499 1.2288603 0.2522424712 0.82 0.4100 ALD1 -0.09222770075 0.9118975 0.1751831913 -0.53 0.6000 CIT 0.00067445016 1.0006747 0.0003391319 1.99 0.0470 DAGE 0.00697116294 1.0069955 0.0042288651 1.65 0.0990 DSEXMale 0.20116657697 1.2228284 0.1653912376 1.22 0.2200

Random effects Group Variable Std Dev Variance
TRANSPLANT_UNIT Intercept 0.00893572849 0.00007984724

confint.coxme(cox2)

                             2.5 %       97.5 %

LDLT1 -1.191361198938 0.3438523103 BMI -0.012570612815 0.0177016505 RBILIRUBIN -0.007251150718 0.0176100095 RINR -0.079641792050 0.2187753732 RCREAT -0.117141593736 0.1774776992 RREN_SUPPre-tx support -2.757614092748 1.3708961453 RAGE 0.008077998027 0.0412246248 WAITLIST_TIME -0.000299635668 0.0004256299 TX_YR -0.045483169544 0.0656563354 HCV1 -0.057798350760 0.5613626207 HCC1 -0.281791222839 0.3818261876 NASH1 -0.298397777346 0.7105721073 ALD1 -0.442594083386 0.2581386819 CIT -0.000003813729 0.0013527141 DAGE -0.001486567298 0.0154288932 DSEXMale -0.129615898176 0.5319490521

#Uni
cox1 <- coxph(Surv(PSURV,PCENS)~ LDLT, data=UK)
forest_model(cox1)

#Just recipient
cox1 <- coxme(formula= Surv(PSURV_years,PCENS)~LDLT+
              #RSEX+
              BMI+  
              RBILIRUBIN+
              RINR+
              RCREAT+
              RREN_SUP+
              RAGE+
              WAITLIST_TIME+
              TX_YR+
              HCV+
              HCC+
              NASH+
              ALD+
              (1|TRANSPLANT_UNIT),
        data=UK)
cox1

Cox mixed-effects model fit by maximum likelihood Data: UK events, n = 1052, 6140 (455 observations deleted due to missingness) Iterations= 12 63 NULL Integrated Fitted Log-likelihood -8606.923 -8525.364 -8517.52

               Chisq    df p    AIC   BIC

Integrated loglik 163.12 14.00 0 135.12 65.70 Penalized loglik 178.81 17.93 0 142.95 54.05

Model: Surv(PSURV_years, PCENS) ~ LDLT + BMI + RBILIRUBIN + RINR + RCREAT + RREN_SUP + RAGE + WAITLIST_TIME + TX_YR + HCV + HCC + NASH + ALD + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z LDLT1 0.12634361188 1.1346720 0.2928621115 0.43 BMI -0.00434497430 0.9956645 0.0065488845 -0.66 RBILIRUBIN 0.01021981629 1.0102722 0.0050377501 2.03 RINR 0.00008861724 1.0000886 0.0057725146 0.02 RCREAT 0.22575233159 1.2532652 0.0681898259 3.31 RREN_SUPPre-tx support 0.27897102663 1.3217690 0.2456649159 1.14 RAGE 0.02092464308 1.0211451 0.0033606007 6.23 WAITLIST_TIME 0.00042653744 1.0004266 0.0001707659 2.50 TX_YR -0.03293101665 0.9676053 0.0123964167 -2.66 HCV1 0.17183808566 1.1874855 0.0807892999 2.13 HCC1 0.41367085416 1.5123593 0.0735249789 5.63 NASH1 0.33846475309 1.4027923 0.1022330965 3.31 ALD1 -0.03022211292 0.9702300 0.0711517571 -0.42 p LDLT1 0.67000000000 BMI 0.51000000000 RBILIRUBIN 0.04200000000 RINR 0.99000000000 RCREAT 0.00093000000 RREN_SUPPre-tx support 0.26000000000 RAGE 0.00000000048 WAITLIST_TIME 0.01200000000 TX_YR 0.00790000000 HCV1 0.03300000000 HCC1 0.00000001800 NASH1 0.00093000000 ALD1 0.67000000000

Random effects Group Variable Std Dev Variance
TRANSPLANT_UNIT Intercept 0.19094314 0.03645928

confint.coxme(cox1)

                            2.5 %        97.5 %

LDLT1 -0.45938061112 0.7120678349 BMI -0.01744274335 0.0087527947 RBILIRUBIN 0.00014431618 0.0202953164 RINR -0.01145641189 0.0116336464 RCREAT 0.08937267981 0.3621319834 RREN_SUPPre-tx support -0.21235880514 0.7703008584 RAGE 0.01420344177 0.0276458444 WAITLIST_TIME 0.00008500563 0.0007680693 TX_YR -0.05772385011 -0.0081381832 HCV1 0.01025948576 0.3334166855 HCC1 0.26662089632 0.5607208120 NASH1 0.13399856012 0.5429309461 ALD1 -0.17252562711 0.1120814013

#Just donor
cox1 <- coxme(formula= Surv(PSURV_years,PCENS)~LDLT+
              CIT+
              DAGE+
              DSEX+
              TX_YR+
              (1|TRANSPLANT_UNIT),
        data=UK)
cox1

Cox mixed-effects model fit by maximum likelihood Data: UK events, n = 1136, 6529 (66 observations deleted due to missingness) Iterations= 5 23 NULL Integrated Fitted Log-likelihood -9342.578 -9320.271 -9312.962

              Chisq   df               p   AIC   BIC

Integrated loglik 44.61 6.00 0.0000000558650 32.61 2.40 Penalized loglik 59.23 9.77 0.0000000040377 39.70 -9.47

Model: Surv(PSURV_years, PCENS) ~ LDLT + CIT + DAGE + DSEX + TX_YR + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z p LDLT1 0.01895274899 1.0191335 0.3010405209 0.06 0.95000000 CIT -0.00008644641 0.9999136 0.0001927383 -0.45 0.65000000 DAGE 0.00965224088 1.0096990 0.0019291147 5.00 0.00000056 DSEXMale 0.15353163218 1.1659447 0.0603082673 2.55 0.01100000 TX_YR -0.02737715959 0.9729942 0.0114388235 -2.39 0.01700000

Random effects Group Variable Std Dev Variance
TRANSPLANT_UNIT Intercept 0.16737583 0.02801467

confint.coxme(cox1)

            2.5 %        97.5 %

LDLT1 -0.583128293 0.6210337907 CIT -0.000471923 0.0002990302 DAGE 0.005794011 0.0135104703 DSEXMale 0.032915098 0.2741481668 TX_YR -0.050254807 -0.0044995126

#donor and recipient
cox2 <- coxme(formula= Surv(PSURV_years,PCENS)~LDLT+
              #RSEX+
              BMI+  
              RBILIRUBIN+
              RINR+
              RCREAT+
              RREN_SUP+
              RAGE+
              WAITLIST_TIME+
              TX_YR+
              HCV+
              HCC+
              NASH+
              ALD+
              CIT+
              DAGE+
              DSEX+
               (1|TRANSPLANT_UNIT),
        data=UK)
cox2

Cox mixed-effects model fit by maximum likelihood Data: UK events, n = 1046, 6085 (510 observations deleted due to missingness) Iterations= 16 83 NULL Integrated Fitted Log-likelihood -8550.296 -8457.482 -8449.321

               Chisq    df p    AIC   BIC

Integrated loglik 185.63 17.00 0 151.63 67.43 Penalized loglik 201.95 21.01 0 159.94 55.90

Model: Surv(PSURV_years, PCENS) ~ LDLT + BMI + RBILIRUBIN + RINR + RCREAT + RREN_SUP + RAGE + WAITLIST_TIME + TX_YR + HCV + HCC + NASH + ALD + CIT + DAGE + DSEX + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z p LDLT1 0.1862613637 1.2047371 0.3152597592 0.59 0.5500000000 BMI -0.0073422704 0.9926846 0.0066779223 -1.10 0.2700000000 RBILIRUBIN 0.0112300085 1.0112933 0.0050507896 2.22 0.0260000000 RINR 0.0029691851 1.0029736 0.0058647373 0.51 0.6100000000 RCREAT 0.2239583521 1.2510189 0.0683259823 3.28 0.0010000000 RREN_SUPPre-tx support 0.2913084441 1.3381773 0.2458794354 1.18 0.2400000000 RAGE 0.0200886458 1.0202918 0.0033910757 5.92 0.0000000031 WAITLIST_TIME 0.0004402304 1.0004403 0.0001705023 2.58 0.0098000000 TX_YR -0.0386085630 0.9621272 0.0125176971 -3.08 0.0020000000 HCV1 0.1797297516 1.1968939 0.0812424548 2.21 0.0270000000 HCC1 0.4106019312 1.5077251 0.0738168431 5.56 0.0000000270 NASH1 0.3269844847 1.3867800 0.1023755104 3.19 0.0014000000 ALD1 -0.0538695516 0.9475557 0.0714679768 -0.75 0.4500000000 CIT -0.0001846231 0.9998154 0.0002026372 -0.91 0.3600000000 DAGE 0.0090706698 1.0091119 0.0020469162 4.43 0.0000094000 DSEXMale 0.1316862255 1.1407503 0.0645943784 2.04 0.0410000000

Random effects Group Variable Std Dev Variance TRANSPLANT_UNIT Intercept 0.2022269 0.0408957

confint.coxme(cox2)

                            2.5 %        97.5 %

LDLT1 -0.44425815480 0.8167808821 BMI -0.02069811493 0.0060135741 RBILIRUBIN 0.00112842931 0.0213315878 RINR -0.00876028942 0.0146986596 RCREAT 0.08730638759 0.3606103166 RREN_SUPPre-tx support -0.20045042664 0.7830673148 RAGE 0.01330649447 0.0268707972 WAITLIST_TIME 0.00009922576 0.0007812350 TX_YR -0.06364395725 -0.0135731688 HCV1 0.01724484203 0.3422146611 HCC1 0.26296824495 0.5582356174 NASH1 0.12223346386 0.5317355055 ALD1 -0.19680550521 0.0890664021 CIT -0.00058989758 0.0002206513 DAGE 0.00497683747 0.0131645022 DSEXMale 0.00249746881 0.2608749823

#bound4 <- bound2 %>% filter(COUNTRY == "US" | COUNTRY == "UK") %>% droplevels()

cox2 <- coxme(formula= Surv(PSURV_years,PCENS)~COUNTRY+
              BMI+  
              RBILIRUBIN+
              RINR+
              RCREAT+
              RAGE+
              WAITLIST_TIME+
              TX_YR+
              CIT+
              DAGE+
              DSEX+
               (1|TRANSPLANT_UNIT),
        data=HCConly)
cox2

Cox mixed-effects model fit by maximum likelihood Data: HCConly events, n = 83, 412 (119 observations deleted due to missingness) Iterations= 10 83 NULL Integrated Fitted Log-likelihood -439.6867 -429.4036 -423.9747

              Chisq    df        p   AIC    BIC

Integrated loglik 20.57 13.00 0.081967 -5.43 -36.88 Penalized loglik 31.42 16.61 0.015243 -1.81 -41.99

Model: Surv(PSURV_years, PCENS) ~ COUNTRY + BMI + RBILIRUBIN + RINR + RCREAT + RAGE + WAITLIST_TIME + TX_YR + CIT + DAGE + DSEX + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z p COUNTRYCAN -1.6247605184 0.1969588 1.0646975105 -1.53 0.1300 COUNTRYUK -0.9772316089 0.3763515 1.0315723031 -0.95 0.3400 BMI 0.0041799280 1.0041887 0.0246082606 0.17 0.8700 RBILIRUBIN 0.0240632878 1.0243551 0.0393129404 0.61 0.5400 RINR -0.2566967389 0.7736028 0.4029856492 -0.64 0.5200 RCREAT 0.4439281176 1.5588184 0.1519691009 2.92 0.0035 RAGE 0.0169760373 1.0171209 0.0149045320 1.14 0.2500 WAITLIST_TIME -0.0001169403 0.9998831 0.0002426874 -0.48 0.6300 TX_YR -0.0281049872 0.9722863 0.0477786972 -0.59 0.5600 CIT -0.0024846183 0.9975185 0.0019311366 -1.29 0.2000 DAGE 0.0231877121 1.0234586 0.0112523078 2.06 0.0390 DSEXMale 0.1126158829 1.1192019 0.2314081284 0.49 0.6300

Random effects Group Variable Std Dev Variance
TRANSPLANT_UNIT Intercept 0.29977763 0.08986663

confint.coxme(cox2)

                  2.5 %       97.5 %

COUNTRYCAN -3.7541555393 0.5046345025 COUNTRYUK -3.0403762150 1.0859129973 BMI -0.0450365933 0.0533964493 RBILIRUBIN -0.0545625929 0.1026891685 RINR -1.0626680373 0.5492745595 RCREAT 0.1399899158 0.7478663194 RAGE -0.0128330267 0.0467851013 WAITLIST_TIME -0.0006023152 0.0003684345 TX_YR -0.1236623817 0.0674524072 CIT -0.0063468915 0.0013776549 DAGE 0.0006830964 0.0456923277 DSEXMale -0.3502003739 0.5754321398

#nonhCC
cox2 <- coxme(formula= Surv(PSURV_years,PCENS)~COUNTRY+
              BMI+  
              RBILIRUBIN+
              RINR+
              RCREAT+
              RAGE+
              WAITLIST_TIME+
              TX_YR+
              CIT+
              DAGE+
              DSEX+
               (1|TRANSPLANT_UNIT),
        data=nonHCC)
cox2

Cox mixed-effects model fit by maximum likelihood Data: nonHCC events, n = 288, 2047 (508 observations deleted due to missingness) Iterations= 12 75 NULL Integrated Fitted Log-likelihood -2037.792 -2009.221 -2003.391

              Chisq    df             p   AIC    BIC

Integrated loglik 57.14 13.00 0.00000016946 31.14 -16.48 Penalized loglik 68.80 17.18 0.00000003947 34.45 -28.46

Model: Surv(PSURV_years, PCENS) ~ COUNTRY + BMI + RBILIRUBIN + RINR + RCREAT + RAGE + WAITLIST_TIME + TX_YR + CIT + DAGE + DSEX + (1 | TRANSPLANT_UNIT) Fixed coefficients coef exp(coef) se(coef) z p COUNTRYCAN -0.66510133495 0.5142214 0.4160242667 -1.60 0.110000 COUNTRYUK 0.30639549620 1.3585195 0.3381680739 0.91 0.360000 BMI -0.04231457737 0.9585682 0.0129733924 -3.26 0.001100 RBILIRUBIN -0.00648612545 0.9935349 0.0137029598 -0.47 0.640000 RINR -0.00816335053 0.9918699 0.0949537812 -0.09 0.930000 RCREAT 0.30125511582 1.3515541 0.0897524503 3.36 0.000790 RAGE 0.02187913711 1.0221202 0.0052581897 4.16 0.000032 WAITLIST_TIME -0.00002206554 0.9999779 0.0001067695 -0.21 0.840000 TX_YR -0.02860797099 0.9717974 0.0221015889 -1.29 0.200000 CIT 0.00016631280 1.0001663 0.0002033148 0.82 0.410000 DAGE 0.01624278227 1.0163754 0.0058167755 2.79 0.005200 DSEXMale 0.06078621735 1.0626717 0.1201231436 0.51 0.610000

Random effects Group Variable Std Dev Variance
TRANSPLANT_UNIT Intercept 0.16078798 0.02585278

confint.coxme(cox2)

                  2.5 %        97.5 %

COUNTRYCAN -1.4971498684 0.1669471985 COUNTRYUK -0.3699406516 0.9827316440 BMI -0.0682613621 -0.0163677926 RBILIRUBIN -0.0338920451 0.0209197942 RINR -0.1980709130 0.1817442119 RCREAT 0.1217502152 0.4807600164 RAGE 0.0113627577 0.0323955165 WAITLIST_TIME -0.0002356044 0.0001914734 TX_YR -0.0728111489 0.0155952069 CIT -0.0002403168 0.0005729424 DAGE 0.0046092313 0.0278763332 DSEXMale -0.1794600699 0.3010325045

#Landmark at 90 days
cox1 <- coxph(formula= Surv(PSURV_years,PCENS)~bound4$COUNTRY+
              bound4$CIT+
              bound4$RSEX+
              bound4$BMI+
              bound4$RBILIRUBIN+
              bound4$RINR+
              bound4$RCREAT+
              bound4$RREN_SUP+
              bound4$DAGE+
              bound4$RAGE+
              bound4$WAITLIST_TIME+
              bound4$DSEX+
              bound4$TX_YR+
              bound4$HCC+
              bound4$PSC+
              bound4$ALD+
              bound4$HCV+
              bound4$Centerpercentile,
              subset = PSURV_years <= 0.2464066,
        data=bound4)
forest_model(cox1)

## Warning in recalculate_width_panels(panel_positions, mapped_text =
## mapped_text, : Unable to resize forest panel to be smaller than its heading;
## consider a smaller text size

#Landmark 90 days to 2 years
cox1 <- coxph(formula= Surv(PSURV_years,PCENS)~bound4$COUNTRY+
              bound4$CIT+
              bound4$RSEX+
              bound4$BMI+
              bound4$RBILIRUBIN+
              bound4$RINR+
              bound4$RCREAT+
              bound4$RREN_SUP+
              bound4$DAGE+
              bound4$RAGE+
              bound4$WAITLIST_TIME+
              bound4$DSEX+
              bound4$TX_YR+
              bound4$HCC+
              bound4$PSC+
              bound4$ALD+
              bound4$HCV+
              bound4$Centerpercentile,
              subset = PSURV_years > 0.2464066 & PSURV_years <= 2,
        data=bound4)
forest_model(cox1)

## Warning in recalculate_width_panels(panel_positions, mapped_text =
## mapped_text, : Unable to resize forest panel to be smaller than its heading;
## consider a smaller text size

#Landmark 5
cox1 <- coxph(formula= Surv(PSURV_years,PCENS)~bound4$COUNTRY+
              bound4$CIT+
              bound4$RSEX+
              bound4$BMI+
              bound4$RBILIRUBIN+
              bound4$RINR+
              bound4$RCREAT+
              bound4$RREN_SUP+
              bound4$DAGE+
              bound4$RAGE+
              bound4$WAITLIST_TIME+
              bound4$DSEX+
              bound4$TX_YR+
              bound4$HCC+
              bound4$PSC+
              bound4$ALD+
              bound4$HCV+
              bound4$Centerpercentile,
              subset = PSURV_years > 2 & PSURV_years <= 5,
        data=bound4)

## Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
## Loglik converged before variable 9 ; coefficient may be infinite.

forest_model(cox1)

## Warning in recalculate_width_panels(panel_positions, mapped_text =
## mapped_text, : Unable to resize forest panel to be smaller than its heading;
## consider a smaller text size

library(scales)
#SEnd Gabi RDS file. 
Trendsovertime <- bound4 %>% select(TX_YR, COUNTRY) 
Trendsovertime <- Trendsovertime %>% group_by(TX_YR, COUNTRY) %>% mutate(count = n())
Trendsovertime <- Trendsovertime %>% group_by(TX_YR) %>% mutate(countyr = n())
Trendsovertime <- Trendsovertime %>% mutate(percentage = count/countyr)

#Counts per year
ggplot(Trendsovertime, aes(factor(TX_YR), y = count, group = COUNTRY, 
                  color = COUNTRY)) +
  geom_line(size = 1.5, alpha = 0.8) +
  geom_point(size = 2) +
  scale_color_npg(name="Country") +
  #scale_color_brewer(name = "Etiology", palette = "Set1")+
  theme_test(base_size = 32) +
  xlab("Year of transplant") +
  ylab("Number of transplants") +
  scale_x_discrete(expand = expansion(mult = c(0, 0))) +
  scale_y_continuous(expand = expansion(mult = c(0, 0)), limits = c(0,350)) +
  annotate("text", x = 8, y = 70, label = "Canada: Cox-Stuart trend test p=0.49", size = 5) +
  annotate("text", x = 8, y = 20, label = "UK: Cox-Stuart trend test p=0.73", size = 5) +
  annotate("text", x = 6, y = 280, label = "US: Cox-Stuart trend test p=0.08", size = 5)

COUNTYEAR <- bound3 %>% group_by(TX_YR, COUNTRY) %>% mutate(countYEAR = n()) %>% ungroup()
COUNTYEAR <- COUNTYEAR %>% select(TX_YR, COUNTRY, countYEAR)

LDLT <- bound5 %>% filter(LDLT=="1" & GRAFT_TYPE == "Segmental") %>% group_by(TX_YR, COUNTRY) %>% mutate(countLDLT = n()) %>% ungroup()
LDLT <- LDLT %>% select(TX_YR, COUNTRY, countLDLT) 
LDLT <- distinct(TX_YR, COUNTRY)

unique_rows1 <- !duplicated(LDLT[c("TX_YR","COUNTRY")])
unique.df1 <- LDLT[unique_rows1,]

unique_rows2 <- !duplicated(COUNTYEAR[c("TX_YR","COUNTRY")])
unique.df2 <- COUNTYEAR[unique_rows2,]

unique.df2
unique.df1
test <- merge(unique.df2, unique.df1, by=c("TX_YR", "COUNTRY"))
test

ggplot(test, aes(factor(TX_YR), y = countLDLT, group = COUNTRY, 
                  color = COUNTRY)) +
  geom_area(size = 1.5, alpha = 0.8, color_palette()) +
  geom_area(data=test, aes(factor(TX_YR), y=countYEAR, group = COUNTRY, color=COUNTRY), size = 1.5, alpha=0.1) +
  scale_color_jama(name="Country") +
  theme_test(base_size = 18) +
  xlab("Year of transplant") +
  ylab("Number of transplants") +
  facet_grid(.~COUNTRY, scales="free_y")+
  scale_x_discrete(expand = expansion(mult = c(0, 0)), breaks=seq(2008, 2018, 4)) +
  scale_y_continuous(expand = expansion(mult = c(0, 0)), limits = c(0, 6660))

#saveRDS(bound3, file = "/Users/Ivanics/Desktop/Research/104. UK vs. US. vs. CAN LDLT/Analysis/bound3forgraph.rds")

#Cox-stuart trend test Canada
library(trend)
Canada <- Trendsovertime %>% filter(COUNTRY == "CAN") %>% select(count, TX_YR)

x <- Canada %>% distinct(TX_YR, .keep_all = T)  %>% arrange((TX_YR)) %>% as.data.frame()
x
x <- x[, "count"]
bartels.test(x)
cs.test(x)

#Canada proportion
x <- c(0.201, 0.172, 0.164, 0.150, 0.185, 0.164, 0.171, 0.156, 0.113, 0.110, 0.138)
cs.test(x)

#Cox-stuart trend test UK
UK <- Trendsovertime %>% filter(COUNTRY == "UK") %>% select(count, TX_YR)

x <- UK %>% distinct(TX_YR, .keep_all = T)  %>% arrange((TX_YR)) %>% as.data.frame()
x
x <- x[, "count"]
bartels.test(x)
cs.test(x)

x <- c(0.0169, 0.0218, 0.0123, 0.0191, 0.0178, 0.0178, 0.0119, 0.024, 0.0148, 0.00824, 0.00389)
cs.test(x)

#Cox-stuart trend test US
US <- Trendsovertime %>% filter(COUNTRY == "US") %>% select(count, TX_YR)

x <- US %>% distinct(TX_YR, .keep_all = T)  %>% arrange((TX_YR)) %>% as.data.frame()
x
x <- x[, "count"]
bartels.test(x)
cs.test(x)

x <- c(0.0341, 0.0324, 0.0419, 0.0366, 0.0381, 0.0406, 0.0418, 0.0491, 0.0451,0.0456, 0.05)
cs.test(x)

Trendsovertimedonortype <- bound3 %>% select(TX_YR, DTYPE, COUNTRY) 

Trendsovertimedonortype <- Trendsovertimedonortype %>% group_by(TX_YR, COUNTRY, DTYPE) %>% mutate(count = n()) %>% ungroup()

Trendsovertimedonortype <- Trendsovertimedonortype %>% group_by(TX_YR, COUNTRY) %>% mutate(countyr = n()) %>% ungroup()

Trendsovertimedonortype <- Trendsovertimedonortype %>% mutate(percentage = count/countyr)

Trendsovertimedonortype <- Trendsovertimedonortype %>% filter(DTYPE == "LDLT")

Trendsovertimedonortype %>% group_by(COUNTRY, TX_YR, percentage) %>% count() %>% ungroup() %>% print(n=40)

Percentage <- c(0.0341, 0.0324, 0.0419, 0.0366, 0.0381, 0.0406, 0.0418, 0.0491, 0.0451,0.0456, 0.05, 0.201, 0.172, 0.164, 0.150, 0.185, 0.164, 0.171, 0.156, 0.113, 0.110, 0.138, 0.0169, 0.0218, 0.0123, 0.0191, 0.0178, 0.0178, 0.0119, 0.024, 0.0148, 0.00824, 0.00389)
Country <- c("US", "US", "US", "US", "US", "US", "US", "US", "US", "US", "US", "CAN", "CAN", "CAN", "CAN", "CAN" , "CAN", "CAN", "CAN", "CAN", "CAN", "CAN", "UK", "UK", "UK", "UK", "UK", "UK", "UK", "UK", "UK", "UK", "UK")
Year <- c(2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018)

Trendsovertimedonortype <- data_frame(Percentage, Country, Year)

Trendsovertimedonortype$Country <- factor(Trendsovertimedonortype$Country)
Trendsovertimedonortype$Country <- relevel(Trendsovertimedonortype$Country, "US")

#Counts per year
ggplot(Trendsovertimedonortype, aes(factor(Year), y = Percentage*100, group = Country, 
                  color = Country)) +
  geom_line(size = 1.5, alpha = 0.8) +
  geom_point(size = 2) +
  scale_color_npg(name="Country") +
  #scale_color_brewer(name = "Etiology", palette = "Set1")+
  theme_test(base_size = 32) +
  xlab("Year of transplant") +
  ylab("Percent of transplants (%)") +
  scale_x_discrete(expand = expansion(mult = c(0, 0))) +
  scale_y_continuous(expand = expansion(mult = c(0, 0)), limits = c(0,25)) +
  annotate("text", x = 8, y = 19, label = "Canada: Cox-Stuart trend test p=0.08", size = 5) +
  annotate("text", x = 8, y = 3, label = "UK: Cox-Stuart trend test p=0.49", size = 5) +
  annotate("text", x = 8, y = 6, label = "US: Cox-Stuart trend test p=0.08", size = 5)

Multiple imputations

library(mice)

UKvUSvCANvsUHNforimputation <- bound2 %>% select(TX_YR, RAGE, PSURV, PCENS, DAGE, DTYPE, DBMI, DCMV, DSEX, BLD_GP_MATCH, GRAFT_TYPE, CIT, RSEX, RETHNIC, BMI, WAITLIST_TIME, TRANSPLANT_UNIT, MELD, RREN_SUP, RVENT, RAB_SURGERY, RLIFE, RASCITES, RENCEPH, RBG, RANTI_HCV, RALBUMIN, RINR, RBILIRUBIN, RCREAT, COUNTRY, UKT_PLDGRP, NASH, PSURV_years, PSURV_1year, PCENS_1year, PSURV_3year, PCENS_3year, PSURV_5year, PCENS_5year, UHN)

sapply(UKvUSvCANvsUHNforimputation, function(x) sum(is.na(x)))

#Imputation
init1 <- mice(UKvUSvCANvsUHNforimputation, m=20, maxit=20, seed = 999)

init1$method

imputed <- complete(init1, 3)
summary(init1)
#densityplot(init1)

sapply(UKvUSvCANvsUHNforimputation, function(x) sum(is.na(x)))

modelFitpre <- with(UKvUSvCANvsUHNforimputation, 
              coxph(formula= Surv(PSURV,PCENS)~COUNTRY+
              CIT+
              RSEX+
              BMI+
              RINR+
              RBILIRUBIN+
              RCREAT+
              RREN_SUP+
              DAGE+
              RAGE+
              WAITLIST_TIME+
              DSEX+
              GRAFT_TYPE+
              TX_YR, data=UKvUSvCANvsUHNforimputation))
summary(modelFitpre)




#Uni
coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY))
summary(pool(coximpute), conf.int = T, exponentiate = T)

coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY, subset =  PSURV_years <= 0.2464066))
summary(pool(coximpute), conf.int = T, exponentiate = T)
                  
coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY, subset =  PSURV_years > 0.2464066 & PSURV_years <= 2))
summary(pool(coximpute), conf.int = T, exponentiate = T)
                                    
coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY, subset =  PSURV_years > 2 & PSURV_years <= 5))
summary(pool(coximpute), conf.int = T, exponentiate = T)
    

           

#Just recipient
coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY+
              DAGE+
              DSEX+
              CIT+
              RSEX+
              BMI+
              RINR+
              RBILIRUBIN+
              RCREAT+
              RREN_SUP+
              RAGE+
              WAITLIST_TIME+
              TX_YR))
summary(pool(coximpute), conf.int = T, exponentiate = T)

coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY+
              DAGE+
              DSEX+
              CIT+
              RSEX+
              BMI+
              RINR+
              RBILIRUBIN+
              RCREAT+
              RREN_SUP+
              RAGE+
              WAITLIST_TIME+
              TX_YR,
              subset =  PSURV_years <= 0.2464066))
summary(pool(coximpute), conf.int = T, exponentiate = T)

coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY+
              DAGE+
              DSEX+
              CIT +
              RSEX+
              BMI+
              RINR+
              RBILIRUBIN+
              RCREAT+
              RREN_SUP+
              RAGE+
              WAITLIST_TIME+
              TX_YR,
              subset =  PSURV_years > 0.2464066 & PSURV_years <= 2))
summary(pool(coximpute), conf.int = T, exponentiate = T)

coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY+
              DAGE+
              DSEX+
              CIT+
              RSEX+
              BMI+
              RINR+
              RBILIRUBIN+
              RCREAT+
              RREN_SUP+
              RAGE+
              WAITLIST_TIME+
              TX_YR,
              subset = PSURV_years > 2 & PSURV_years <= 5))
summary(pool(coximpute), conf.int = T, exponentiate = T)



#Just donor
coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY+
              CIT+
              DAGE+
              DSEX+
              TX_YR))
summary(pool(coximpute), conf.int = T, exponentiate = T)

coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY+
              CIT+
              DAGE+
              DSEX+
              TX_YR,
              subset =  PSURV_years <= 0.2464066))
summary(pool(coximpute), conf.int = T, exponentiate = T)

coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY+
              CIT+
              DAGE+
              DSEX+
              TX_YR,
              subset =  PSURV_years > 0.2464066 & PSURV_years <= 2))
summary(pool(coximpute), conf.int = T, exponentiate = T)

coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY+
              CIT+
              DAGE+
              DSEX+
              TX_YR, 
              subset = PSURV_years > 2 & PSURV_years <= 5))
summary(pool(coximpute), conf.int = T, exponentiate = T)




#Overall
coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY+
              CIT+
              RSEX+
              BMI+
              RINR+
              RBILIRUBIN+
              RCREAT+
              RREN_SUP+
              DAGE+
              RAGE+
              WAITLIST_TIME+
              DSEX+
              TX_YR))
summary(pool(coximpute), conf.int = T, exponentiate = T)
  
#Landmark 0-90 days
coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY+
              CIT+
              RSEX+
              BMI+
              RINR+
              RBILIRUBIN+
              RCREAT+
              RREN_SUP+
              DAGE+
              RAGE+
              WAITLIST_TIME+
              DSEX+
              TX_YR, 
              subset =  PSURV_years <= 0.2464066))
summary(pool(coximpute), conf.int = T, exponentiate = T)

#Landmark 90-days to 2 years
coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY+
              CIT+
              RSEX+
              BMI+
              RINR+
              RBILIRUBIN+
              RCREAT+
              RREN_SUP+
              DAGE+
              RAGE+
              WAITLIST_TIME+
              DSEX+
              TX_YR, 
              subset = PSURV_years > 0.2464066 & PSURV_years <= 2))
summary(pool(coximpute), conf.int = T, exponentiate = T)

#Landmark 2years to 5 years
coximpute <- with(init1, coxph(Surv(PSURV,PCENS)~ COUNTRY+
              CIT+
              RSEX+
              BMI+
              RINR+
              RBILIRUBIN+
              RCREAT+
              RREN_SUP+
              DAGE+
              RAGE+
              WAITLIST_TIME+
              DSEX+
              TX_YR, 
              subset = PSURV_years > 2 & PSURV_years <= 5))
summary(pool(coximpute), conf.int = T, exponentiate = T)